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.
###############################################################################
Custom Compressor Unit Model#
Author: Andrew Lee
Maintainer: Andrew Lee
Updated: 2023-06-01
To demonstrate creation of a new unit model, we will create a constant-heat-capacity ideal-gas isentropic compressor. This will be a simple textbook model. We will utilize the mass and energy balances provided by IDAES control volumes, but we will write our own isentropic constraint based off of equations 7.18 and 7.23 from “Introduction to Chemical Engineering Thermodynamics” by J.M. Smith, H.C. Van Ness, and M.M. Abbott.
The outlet temperature of an ideal gas undergoing isentropic compression is given by $\( t_{out} = t_{in} + \frac{1}{\eta} \left(t_{in} \left(\frac{p_{out}}{p_{in}}\right)^{\frac{\gamma - 1}{\gamma}} - t_{in}\right) \)\( where \)p\( is pressure, \)t\( is temperature, and \)\gamma$ is the ratio of constant pressure heat capacity to constant volume heat capacity.
We will begin with relevant imports. We will need
Pyomo for writing our energy balance constraints
ConfigBlocks for specifying options for our compressor
ControlVolume0DBlocks for creating the appropriate state blocks for the inlet and outlet and for defining mas balances
IdealParameterBlock which provides a simple ideal-gas property package.
A few other helpful functions and enums from IDAES
import pyomo.environ as pe
from pyomo.common.config import ConfigBlock, ConfigValue, In
from idaes.core import (
ControlVolume0DBlock,
declare_process_block_class,
EnergyBalanceType,
MomentumBalanceType,
MaterialBalanceType,
UnitModelBlockData,
useDefault,
FlowsheetBlock,
)
from idaes.core.util.config import is_physical_parameter_block
from idaes_examples.mod.methanol.methanol_param_VLE import PhysicalParameterBlock
from idaes.core.util.misc import add_object_reference
Now, we can write a function to create a control volume for our compressor. The control volume will define the inlet and outlet streams along with the appropriate state variables (specified by the property package). We will also use the control volume to create mass and energy balance constraints.
Our function will take the compressor unit model object, the name of the control volume, and configuration options as arguments. Our compressor will only support steady-state models, so we will first ensure that dynamic and has_holdup are both False.
Next, we will create a 0D control volume. We are using a 0D control volume because our model does not depend on space. We then
Attach the control volume to the compressor
Create the appropriate state blocks with the control volume (for the inlet and outlet streams)
Use the control volume to add mass balance constraints
Use the control volume to add energy balance constraints
def make_control_volume(unit, name, config):
if config.dynamic is not False:
raise ValueError("IdealGasIsentropcCompressor does not support dynamics")
if config.has_holdup is not False:
raise ValueError("IdealGasIsentropcCompressor does not support holdup")
control_volume = ControlVolume0DBlock(
property_package=config.property_package,
property_package_args=config.property_package_args,
)
setattr(unit, name, control_volume)
control_volume.add_state_blocks(has_phase_equilibrium=config.has_phase_equilibrium)
control_volume.add_material_balances(
balance_type=config.material_balance_type,
has_phase_equilibrium=config.has_phase_equilibrium,
)
control_volume.add_total_enthalpy_balances(
has_heat_of_reaction=False, has_heat_transfer=False, has_work_transfer=True
)
Next, we will write a function to add constraints to specify that the compressor is isentropic.
Create a
pressure_ratiovariable to represent \(p_{out}/p_{in}\). The lower bound is \(1\), because we only want to allow compression (and not expansion).Create a
ConstraintListto hold the constraints.Add the
ConstraintListto the compressorCreate the local variables
inletandoutletto reference the inlet and outlet state blocks.Add a constraint relating the inlet pressure, outlet pressure, and pressure ratio variables: \begin{align} p_{in} p_{ratio} = p_{out} \end{align}
Add a constraint relating the inlet and outlet temperatures: \begin{align} & t_{out} = t_{in} + \frac{1}{\eta} \left(t_{in} p_{ratio}^{\frac{\gamma - 1}{\gamma}} - t_{in}\right) \end{align}
def add_isentropic(unit, name, config):
unit.pressure_ratio = pe.Var(initialize=1.0, bounds=(1, None))
cons = pe.ConstraintList()
setattr(unit, name, cons)
inlet = unit.control_volume.properties_in[0.0]
outlet = unit.control_volume.properties_out[0.0]
gamma = inlet.params.gamma
cons.add(inlet.pressure * unit.pressure_ratio == outlet.pressure)
cons.add(
outlet.temperature
== (
inlet.temperature
+ 1
/ config.compressor_efficiency
* (
inlet.temperature * unit.pressure_ratio ** ((gamma - 1) / gamma)
- inlet.temperature
)
)
)
We also need a function to specify configuration options for the compressor.
def make_compressor_config_block(config):
config.declare(
"material_balance_type",
ConfigValue(
default=MaterialBalanceType.componentPhase, domain=In(MaterialBalanceType)
),
)
config.declare(
"energy_balance_type",
ConfigValue(
default=EnergyBalanceType.enthalpyTotal,
domain=In([EnergyBalanceType.enthalpyTotal]),
),
)
config.declare(
"momentum_balance_type",
ConfigValue(
default=MomentumBalanceType.none, domain=In([MomentumBalanceType.none])
),
)
config.declare(
"has_phase_equilibrium", ConfigValue(default=False, domain=In([False]))
)
config.declare(
"has_pressure_change", ConfigValue(default=False, domain=In([False]))
)
config.declare(
"property_package",
ConfigValue(default=useDefault, domain=is_physical_parameter_block),
)
config.declare("property_package_args", ConfigBlock(implicit=True))
config.declare("compressor_efficiency", ConfigValue(default=0.75, domain=float))
Finally, we can define the ideal-gas isentropic compressor. To do so, we create a class called IdealGasIsentropicCompressorData and use the declare_process_block_class decorator. For now, just consider the decorator to be boiler-plate. We then need to define the config block and write the build method. The build method should always call super. Next, we simply call the functions we wrote to build the control volume, energy balance, and electricity requirement performance equation. Finally, we need to call self.add_inlet_port() and self.add_outlet_port(). These methods need to be called in order to create the ports which are used for connecting the unit to other units.
@declare_process_block_class("IdealGasIsentropicCompressor")
class IdealGasIsentropicCompressorData(UnitModelBlockData):
CONFIG = UnitModelBlockData.CONFIG()
make_compressor_config_block(CONFIG)
def build(self):
super(IdealGasIsentropicCompressorData, self).build()
make_control_volume(self, "control_volume", self.config)
add_isentropic(self, "isentropic", self.config)
self.add_inlet_port()
self.add_outlet_port()
add_object_reference(self, "work", self.control_volume.work[0.0])
The compressor model is complete and can now be used like other IDAES unit models. Note that the input temperature is in hectoKelvin, the input pressure is in MPa and energy units are in MJ. This is to simplify user input and is accounted for in the property package files; the standard unit definitions may be found in the metadata section at the end of the main parameter property package.
m = pe.ConcreteModel()
m.fs = FlowsheetBlock(dynamic=False)
m.fs.properties = props = PhysicalParameterBlock(
Cp=0.038056, valid_phase="Vap"
) # MJ/kmol-K
m.fs.compressor = IdealGasIsentropicCompressor(
property_package=props, has_phase_equilibrium=False
)
m.fs.compressor.inlet.flow_mol.fix(1) # kmol
m.fs.compressor.inlet.mole_frac_comp[0, "CH3OH"].fix(0.25)
m.fs.compressor.inlet.mole_frac_comp[0, "CH4"].fix(0.25)
m.fs.compressor.inlet.mole_frac_comp[0, "H2"].fix(0.25)
m.fs.compressor.inlet.mole_frac_comp[0, "CO"].fix(0.25)
m.fs.compressor.inlet.pressure.fix(0.14) # MPa
m.fs.compressor.inlet.temperature.fix(2.9315) # hK [100K]
m.fs.compressor.outlet.pressure.fix(0.56) # MPa
opt = pe.SolverFactory("ipopt")
opt.options["linear_solver"] = "ma27"
res = opt.solve(m, tee=True)
print(res.solver.termination_condition)
m.fs.compressor.outlet.display()
print("work: ", round(m.fs.compressor.work.value, 2), " MJ") # MJ
Ipopt 3.13.2: linear_solver=ma27
******************************************************************************
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...: 54
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 14
Total number of variables............................: 18
variables with only lower bounds: 5
variables with lower and upper bounds: 12
variables with only upper bounds: 0
Total number of equality constraints.................: 18
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 5.00e-01 1.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 0.0000000e+00 4.71e+00 2.89e+02 -1.0 4.97e+00 - 3.41e-03 1.00e+00f 1
2 0.0000000e+00 2.66e-15 1.00e-06 -1.0 1.16e+00 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 2
(scaled) (unscaled)
Objective...............: 0.0000000000000000e+00 0.0000000000000000e+00
Dual infeasibility......: 0.0000000000000000e+00 0.0000000000000000e+00
Constraint violation....: 2.6645352591003757e-15 2.6645352591003757e-15
Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00
Overall NLP error.......: 2.6645352591003757e-15 2.6645352591003757e-15
Number of objective function evaluations = 3
Number of objective gradient evaluations = 3
Number of equality constraint evaluations = 3
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 3
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 2
Total CPU secs in IPOPT (w/o function evaluations) = 0.001
Total CPU secs in NLP function evaluations = 0.000
EXIT: Optimal Solution Found.
optimal
outlet : Size=1
Key : Name : Value
None : flow_mol : {0.0: 1.0}
: mole_frac_comp : {(0.0, 'CH3OH'): 0.25, (0.0, 'CH4'): 0.25, (0.0, 'CO'): 0.25, (0.0, 'H2'): 0.25}
: pressure : {0.0: 0.56}
: temperature : {0.0: 4.314183563052119}
work: 5.26 MJ