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.
###############################################################################
Flowsheet Gibbs Reactor Simulation and Optimization of Steam Methane Reforming#
Author: Brandon Paul
Maintainer: Brandon Paul
Updated: 2023-06-01
Learning Outcomes#
Call and implement the IDAES GibbsReactor unit model
Construct a steady-state flowsheet using the IDAES unit model library
Connecting unit models in a flowsheet using Arcs
Fomulate and solve an optimization problem
Defining an objective function
Setting variable bounds
Adding additional constraints
Problem Statement#
Following the previous example of Steam Methane Reformation in an Equilibrium Reactor, this example solves the flowsheet using a Gibbs Reactor instead. The steam methane reformation example is adapted from S.Z. Abbas, V. Dupont, T. Mahmud, Kinetics study and modelling of steam methane reforming process over a NiO/Al2O3 catalyst in an adiabatic packed bed reactor. Int. J. Hydrogen Energy, 42 (2017), pp. 2889-2903. Typically, the process follows the chemical equations below:
CH4 + H2O → CO + 3H2
CO + H2O → CO2 + H2
However, the GibbsReactor unit model solves the equilibrium by minimizing Gibbs free energy. Conveniently, this eliminates the need for a reaction package although a thermophysical package is still required.
The flowsheet that we will be using for this module is shown below with the stream conditions. As in the prior example, we will be processing natural gas and steam feeds of fixed composition to produce hydrogen. The process consists of a mixer M101 for the two inlet streams, a compressor to compress the feed to the reaction pressure, a heater H101 to heat the feed to the reaction temperature, and a GibbsReactor unit R101. We will use thermophysical properties following the Peng-Robinsion cubic equation of state for this flowsheet.
The state variables chosen for the property package are total molar flows of each stream, temperature of each stream and pressure of each stream, and mole fractions of each component in each stream. The components considered are: CH4, H2O, CO, CO2, and H2 and the process occurs in vapor phase only. Therefore, every stream has 1 flow variable, 5 mole fraction variables, 1 temperature and 1 pressure variable.
Importing Required Pyomo and IDAES Components#
To construct a flowsheet, we will need several components from the Pyomo and IDAES packages, as well as some utility tools to build the flowsheet. For further details on these components, please refer to the Pyomo documentation.
From IDAES, we will be needing the FlowsheetBlock and the following unit models:
Feed
Mixer
Compressor
Heater
GibbsReactor
Product
from pyomo.environ import (
Constraint,
Var,
ConcreteModel,
Expression,
Objective,
TransformationFactory,
value,
units as pyunits,
)
from pyomo.network import Arc
from idaes.core import FlowsheetBlock
from idaes.models.properties.modular_properties import GenericParameterBlock
from idaes.models.unit_models import (
Feed,
Mixer,
Compressor,
Heater,
GibbsReactor,
Product,
)
from idaes.core.solvers import get_solver
from idaes.core.util.model_statistics import degrees_of_freedom
from idaes.core.util.initialization import propagate_state
Importing Required Thermophysical Package#
As mentioned earlier, the GibbsReactor does not require a reaction package.
Let us import the following module from the IDAES library:
natural_gas_PR as get_prop (method to get configuration dictionary)
from idaes.models_extra.power_generation.properties.natural_gas_PR import get_prop
Constructing the Flowsheet#
Let us create a ConcreteModel and add the flowsheet block.
m = ConcreteModel()
m.fs = FlowsheetBlock(dynamic=False)
We now need to add the property packages to the flowsheet. Unlike the previous example, we do not need to add a reaction package for the reactor model to calculate results. We will use the Modular Property Framework to build a state block for the parameter dictionary.
thermo_props_config_dict = get_prop(components=["CH4", "H2O", "H2", "CO", "CO2"])
m.fs.thermo_params = GenericParameterBlock(**thermo_props_config_dict)
Adding Unit Models#
Let us start adding the unit models we have imported to the flowsheet:
m.fs.CH4 = Feed(property_package=m.fs.thermo_params)
m.fs.H2O = Feed(property_package=m.fs.thermo_params)
m.fs.PROD = Product(property_package=m.fs.thermo_params)
m.fs.M101 = Mixer(
property_package=m.fs.thermo_params, inlet_list=["methane_feed", "steam_feed"]
)
m.fs.H101 = Heater(
property_package=m.fs.thermo_params,
has_pressure_change=False,
has_phase_equilibrium=False,
)
m.fs.C101 = Compressor(property_package=m.fs.thermo_params)
ERROR: Rule failed for Expression
'fs.M101.methane_feed_state[0.0].compress_fact_phase' with index Vap:
RuntimeError: Cubic root external functions are not available.
ERROR: Constructing component
'fs.M101.methane_feed_state[0.0].compress_fact_phase' from data=None failed:
RuntimeError: Cubic root external functions are not available.
ERROR: Rule failed for Expression
'fs.M101.methane_feed_state[0.0].enth_mol_phase' with index Vap: RuntimeError:
Cubic root external functions are not available.
ERROR: Constructing component 'fs.M101.methane_feed_state[0.0].enth_mol_phase'
from data=None failed:
RuntimeError: Cubic root external functions are not available.
ERROR: Rule failed for Expression
'fs.M101.methane_feed_state[0.0]._enthalpy_flow_term' with index Vap:
RuntimeError: Cubic root external functions are not available.
ERROR: Constructing component
'fs.M101.methane_feed_state[0.0]._enthalpy_flow_term' from data=None failed:
RuntimeError: Cubic root external functions are not available.
ERROR: Rule failed when generating expression for Constraint
fs.M101.enthalpy_mixing_equations with index 0.0: RuntimeError: Cubic root
external functions are not available.
ERROR: Constructing component 'fs.M101.enthalpy_mixing_equations' from
data=None failed:
RuntimeError: Cubic root external functions are not available.
2024-05-09 18:40:09 [ERROR] idaes.core.base.process_block: Failure in build: fs.M101
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/property_base.py", line 792, in __getattr__
return super().__getattr__(attr)
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/core/base/block.py", line 550, in __getattr__
raise AttributeError(
AttributeError: 'GenericStateBlockData' object has no attribute '_enthalpy_flow_term'
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/models/properties/modular_properties/state_definitions/FTPx.py", line 279, in get_enthalpy_flow_terms_FTPx
eflow = b._enthalpy_flow_term
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/core/base/property_base.py", line 794, in __getattr__
return build_on_demand(self, attr)
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/core/base/util.py", line 85, in build_on_demand
raise PropertyPackageError(
idaes.core.util.exceptions.PropertyPackageError: fs.M101.methane_feed_state[0.0] _enthalpy_flow_term does not exist, but is a protected attribute. Check the naming of your components to avoid any reserved names
During handling of the above exception, another exception occurred:
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/property_base.py", line 792, in __getattr__
return super().__getattr__(attr)
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/core/base/block.py", line 550, in __getattr__
raise AttributeError(
AttributeError: 'GenericStateBlockData' object has no attribute 'enth_mol_phase'
During handling of the above exception, another exception occurred:
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/property_base.py", line 792, in __getattr__
return super().__getattr__(attr)
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/core/base/block.py", line 550, in __getattr__
raise AttributeError(
AttributeError: 'GenericStateBlockData' object has no attribute 'compress_fact_phase'
During handling of the above exception, another exception occurred:
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/unit_models/mixer.py", line 450, in build
self.add_energy_mixing_equations(
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/models/unit_models/mixer.py", line 797, in add_energy_mixing_equations
def enthalpy_mixing_equations(b, t):
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/core/base/block.py", line 77, in __call__
setattr(
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/core/base/block.py", line 571, in __setattr__
self.add_component(name, val)
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/core/base/block.py", line 1129, in add_component
val.construct(data)
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/core/base/constraint.py", line 720, in construct
self._setitem_when_not_present(index, rule(block, index))
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/core/base/initializer.py", line 316, in __call__
return self._fcn(parent, idx)
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/models/unit_models/mixer.py", line 799, in enthalpy_mixing_equations
sum(
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/models/unit_models/mixer.py", line 800, in <genexpr>
sum(
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/models/unit_models/mixer.py", line 801, in <genexpr>
inlet_blocks[i][t].get_enthalpy_flow_terms(p)
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/models/properties/modular_properties/state_definitions/FTPx.py", line 285, in get_enthalpy_flow_terms_FTPx
eflow = b._enthalpy_flow_term = Expression(b.phase_list, rule=rule_eflow)
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/core/base/block.py", line 571, in __setattr__
self.add_component(name, val)
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/core/base/block.py", line 1129, in add_component
val.construct(data)
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/core/base/expression.py", line 375, in construct
self._construct_from_rule_using_setitem()
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/core/base/indexed_component.py", line 784, in _construct_from_rule_using_setitem
self._setitem_when_not_present(index, rule(block, index))
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/core/base/initializer.py", line 316, in __call__
return self._fcn(parent, idx)
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/models/properties/modular_properties/state_definitions/FTPx.py", line 283, in rule_eflow
return b.flow_mol_phase[p] * b.enth_mol_phase[p]
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/core/base/property_base.py", line 794, in __getattr__
return build_on_demand(self, attr)
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/core/base/util.py", line 220, in build_on_demand
f()
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/models/properties/modular_properties/base/generic_property.py", line 3537, in _enth_mol_phase
self.enth_mol_phase = Expression(self.phase_list, rule=rule_enth_mol_phase)
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/core/base/block.py", line 571, in __setattr__
self.add_component(name, val)
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/core/base/block.py", line 1129, in add_component
val.construct(data)
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/core/base/expression.py", line 375, in construct
self._construct_from_rule_using_setitem()
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/core/base/indexed_component.py", line 784, in _construct_from_rule_using_setitem
self._setitem_when_not_present(index, rule(block, index))
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/core/base/initializer.py", line 316, in __call__
return self._fcn(parent, idx)
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/models/properties/modular_properties/base/generic_property.py", line 3535, in rule_enth_mol_phase
return p_config.equation_of_state.enth_mol_phase(b, p)
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/models/properties/modular_properties/eos/ceos.py", line 690, in enth_mol_phase
Z = blk.compress_fact_phase[p]
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/core/base/property_base.py", line 794, in __getattr__
return build_on_demand(self, attr)
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/core/base/util.py", line 220, in build_on_demand
f()
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/models/properties/modular_properties/base/generic_property.py", line 3164, in _compress_fact_phase
self.compress_fact_phase = Expression(
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/core/base/block.py", line 571, in __setattr__
self.add_component(name, val)
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/core/base/block.py", line 1129, in add_component
val.construct(data)
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/core/base/expression.py", line 375, in construct
self._construct_from_rule_using_setitem()
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/core/base/indexed_component.py", line 784, in _construct_from_rule_using_setitem
self._setitem_when_not_present(index, rule(block, index))
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/core/base/initializer.py", line 316, in __call__
return self._fcn(parent, idx)
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/models/properties/modular_properties/base/generic_property.py", line 3162, in rule_Z_phase
return p_config.equation_of_state.compress_fact_phase(b, p)
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/models/properties/modular_properties/eos/ceos.py", line 551, in compress_fact_phase
expr_write = CubicThermoExpressions(b)
File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/models/properties/modular_properties/eos/ceos_common.py", line 106, in __init__
raise RuntimeError("Cubic root external functions are not available.")
RuntimeError: Cubic root external functions are not available.
ERROR: Constructing component 'fs.M101' from data=None failed:
RuntimeError: Cubic root external functions are not available.
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/core/base/property_base.py:792, in StateBlockData.__getattr__(self, attr)
788 try:
789 # try the Pyomo Block's __getattr__ method first which will return
790 # decorators for creating components on the block (e.g. Expression,
791 # Constraint, ...).
--> 792 return super().__getattr__(attr)
793 except AttributeError:
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/core/base/block.py:550, in BlockData.__getattr__(self, val)
548 # Since the base classes don't support getattr, we can just
549 # throw the "normal" AttributeError
--> 550 raise AttributeError(
551 "'%s' object has no attribute '%s'" % (self.__class__.__name__, val)
552 )
AttributeError: 'GenericStateBlockData' object has no attribute '_enthalpy_flow_term'
During handling of the above exception, another exception occurred:
PropertyPackageError Traceback (most recent call last)
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/models/properties/modular_properties/state_definitions/FTPx.py:279, in define_state.<locals>.get_enthalpy_flow_terms_FTPx(b, p)
278 try:
--> 279 eflow = b._enthalpy_flow_term
280 except AttributeError:
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/core/base/property_base.py:794, in StateBlockData.__getattr__(self, attr)
793 except AttributeError:
--> 794 return build_on_demand(self, attr)
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/core/base/util.py:85, in build_on_demand(self, attr)
82 if attr == "domain" or attr.startswith("_"):
83 # Don't interfere with anything by getting attributes that are
84 # none of my business
---> 85 raise PropertyPackageError(
86 "{} {} does not exist, but is a protected "
87 "attribute. Check the naming of your "
88 "components to avoid any reserved names".format(self.name, attr)
89 )
91 if attr == "config":
PropertyPackageError: fs.M101.methane_feed_state[0.0] _enthalpy_flow_term does not exist, but is a protected attribute. Check the naming of your components to avoid any reserved names
During handling of the above exception, another exception occurred:
AttributeError Traceback (most recent call last)
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/core/base/property_base.py:792, in StateBlockData.__getattr__(self, attr)
788 try:
789 # try the Pyomo Block's __getattr__ method first which will return
790 # decorators for creating components on the block (e.g. Expression,
791 # Constraint, ...).
--> 792 return super().__getattr__(attr)
793 except AttributeError:
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/core/base/block.py:550, in BlockData.__getattr__(self, val)
548 # Since the base classes don't support getattr, we can just
549 # throw the "normal" AttributeError
--> 550 raise AttributeError(
551 "'%s' object has no attribute '%s'" % (self.__class__.__name__, val)
552 )
AttributeError: 'GenericStateBlockData' object has no attribute 'enth_mol_phase'
During handling of the above exception, another exception occurred:
AttributeError Traceback (most recent call last)
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/core/base/property_base.py:792, in StateBlockData.__getattr__(self, attr)
788 try:
789 # try the Pyomo Block's __getattr__ method first which will return
790 # decorators for creating components on the block (e.g. Expression,
791 # Constraint, ...).
--> 792 return super().__getattr__(attr)
793 except AttributeError:
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/core/base/block.py:550, in BlockData.__getattr__(self, val)
548 # Since the base classes don't support getattr, we can just
549 # throw the "normal" AttributeError
--> 550 raise AttributeError(
551 "'%s' object has no attribute '%s'" % (self.__class__.__name__, val)
552 )
AttributeError: 'GenericStateBlockData' object has no attribute 'compress_fact_phase'
During handling of the above exception, another exception occurred:
RuntimeError Traceback (most recent call last)
Cell In[6], line 4
2 m.fs.H2O = Feed(property_package=m.fs.thermo_params)
3 m.fs.PROD = Product(property_package=m.fs.thermo_params)
----> 4 m.fs.M101 = Mixer(
5 property_package=m.fs.thermo_params, inlet_list=["methane_feed", "steam_feed"]
6 )
7 m.fs.H101 = Heater(
8 property_package=m.fs.thermo_params,
9 has_pressure_change=False,
10 has_phase_equilibrium=False,
11 )
12 m.fs.C101 = Compressor(property_package=m.fs.thermo_params)
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/unit_models/mixer.py:450, in MixerData.build(self)
442 raise BurntToast(
443 "{} received unrecognised value for "
444 "material_mixing_type argument. This "
445 "should not occur, so please contact "
446 "the IDAES developers with this bug.".format(self.name)
447 )
449 if self.config.energy_mixing_type == MixingType.extensive:
--> 450 self.add_energy_mixing_equations(
451 inlet_blocks=inlet_blocks, mixed_block=mixed_block
452 )
453 elif self.config.energy_mixing_type == MixingType.none:
454 pass
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/models/unit_models/mixer.py:797, in MixerData.add_energy_mixing_equations(self, inlet_blocks, mixed_block)
791 def add_energy_mixing_equations(self, inlet_blocks, mixed_block):
792 """
793 Add energy mixing equations (total enthalpy balance).
794 """
796 @self.Constraint(self.flowsheet().time, doc="Energy balances")
--> 797 def enthalpy_mixing_equations(b, t):
798 return 0 == (
799 sum(
800 sum(
(...)
809 )
810 )
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/core/base/block.py:77, in _generic_component_decorator.__call__(self, rule)
76 def __call__(self, rule):
---> 77 setattr(
78 self._block,
79 rule.__name__,
80 self._component(*self._args, rule=rule, **(self._kwds)),
81 )
82 return rule
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/constraint.py:720, in Constraint.construct(self, data)
717 else:
718 # Bypass the index validation and create the member directly
719 for index in self.index_set():
--> 720 self._setitem_when_not_present(index, rule(block, index))
721 except Exception:
722 err = sys.exc_info()[1]
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/models/unit_models/mixer.py:799, in MixerData.add_energy_mixing_equations.<locals>.enthalpy_mixing_equations(b, t)
796 @self.Constraint(self.flowsheet().time, doc="Energy balances")
797 def enthalpy_mixing_equations(b, t):
798 return 0 == (
--> 799 sum(
800 sum(
801 inlet_blocks[i][t].get_enthalpy_flow_terms(p)
802 for p in mixed_block.phase_list
803 )
804 for i in range(len(inlet_blocks))
805 )
806 - sum(
807 mixed_block[t].get_enthalpy_flow_terms(p)
808 for p in mixed_block.phase_list
809 )
810 )
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/models/unit_models/mixer.py:800, in <genexpr>(.0)
796 @self.Constraint(self.flowsheet().time, doc="Energy balances")
797 def enthalpy_mixing_equations(b, t):
798 return 0 == (
799 sum(
--> 800 sum(
801 inlet_blocks[i][t].get_enthalpy_flow_terms(p)
802 for p in mixed_block.phase_list
803 )
804 for i in range(len(inlet_blocks))
805 )
806 - sum(
807 mixed_block[t].get_enthalpy_flow_terms(p)
808 for p in mixed_block.phase_list
809 )
810 )
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/models/unit_models/mixer.py:801, in <genexpr>(.0)
796 @self.Constraint(self.flowsheet().time, doc="Energy balances")
797 def enthalpy_mixing_equations(b, t):
798 return 0 == (
799 sum(
800 sum(
--> 801 inlet_blocks[i][t].get_enthalpy_flow_terms(p)
802 for p in mixed_block.phase_list
803 )
804 for i in range(len(inlet_blocks))
805 )
806 - sum(
807 mixed_block[t].get_enthalpy_flow_terms(p)
808 for p in mixed_block.phase_list
809 )
810 )
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/models/properties/modular_properties/state_definitions/FTPx.py:285, in define_state.<locals>.get_enthalpy_flow_terms_FTPx(b, p)
282 def rule_eflow(b, p):
283 return b.flow_mol_phase[p] * b.enth_mol_phase[p]
--> 285 eflow = b._enthalpy_flow_term = Expression(b.phase_list, rule=rule_eflow)
286 return eflow[p]
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/expression.py:375, in Expression.construct(self, data)
372 try:
373 # We do not (currently) accept data for constructing Constraints
374 assert data is None
--> 375 self._construct_from_rule_using_setitem()
376 finally:
377 timer.report()
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/core/base/indexed_component.py:784, in IndexedComponent._construct_from_rule_using_setitem(self)
782 else:
783 for index in self.index_set():
--> 784 self._setitem_when_not_present(index, rule(block, index))
785 except:
786 err = sys.exc_info()[1]
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/models/properties/modular_properties/state_definitions/FTPx.py:283, in define_state.<locals>.get_enthalpy_flow_terms_FTPx.<locals>.rule_eflow(b, p)
282 def rule_eflow(b, p):
--> 283 return b.flow_mol_phase[p] * b.enth_mol_phase[p]
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/core/base/property_base.py:794, in StateBlockData.__getattr__(self, attr)
792 return super().__getattr__(attr)
793 except AttributeError:
--> 794 return build_on_demand(self, attr)
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/core/base/util.py:220, in build_on_demand(self, attr)
218 if callable(f):
219 try:
--> 220 f()
221 except Exception:
222 # Clear call list and reraise error
223 clear_call_list(self, attr)
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/models/properties/modular_properties/base/generic_property.py:3537, in GenericStateBlockData._enth_mol_phase(self)
3534 p_config = b.params.get_phase(p).config
3535 return p_config.equation_of_state.enth_mol_phase(b, p)
-> 3537 self.enth_mol_phase = Expression(self.phase_list, rule=rule_enth_mol_phase)
3538 except AttributeError:
3539 self.del_component(self.enth_mol_phase)
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/expression.py:375, in Expression.construct(self, data)
372 try:
373 # We do not (currently) accept data for constructing Constraints
374 assert data is None
--> 375 self._construct_from_rule_using_setitem()
376 finally:
377 timer.report()
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/core/base/indexed_component.py:784, in IndexedComponent._construct_from_rule_using_setitem(self)
782 else:
783 for index in self.index_set():
--> 784 self._setitem_when_not_present(index, rule(block, index))
785 except:
786 err = sys.exc_info()[1]
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/models/properties/modular_properties/base/generic_property.py:3535, in GenericStateBlockData._enth_mol_phase.<locals>.rule_enth_mol_phase(b, p)
3533 def rule_enth_mol_phase(b, p):
3534 p_config = b.params.get_phase(p).config
-> 3535 return p_config.equation_of_state.enth_mol_phase(b, p)
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/models/properties/modular_properties/eos/ceos.py:690, in Cubic.enth_mol_phase(blk, p)
688 B = getattr(blk, cname + "_B")[p]
689 dam_dT = getattr(blk, cname + "_dam_dT")[p]
--> 690 Z = blk.compress_fact_phase[p]
691 R = Cubic.gas_constant(blk)
692 T = blk.temperature
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/core/base/property_base.py:794, in StateBlockData.__getattr__(self, attr)
792 return super().__getattr__(attr)
793 except AttributeError:
--> 794 return build_on_demand(self, attr)
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/core/base/util.py:220, in build_on_demand(self, attr)
218 if callable(f):
219 try:
--> 220 f()
221 except Exception:
222 # Clear call list and reraise error
223 clear_call_list(self, attr)
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/models/properties/modular_properties/base/generic_property.py:3164, in GenericStateBlockData._compress_fact_phase(self)
3161 p_config = b.params.get_phase(p).config
3162 return p_config.equation_of_state.compress_fact_phase(b, p)
-> 3164 self.compress_fact_phase = Expression(
3165 self.phase_list, doc="Compressibility of each phase", rule=rule_Z_phase
3166 )
3167 except AttributeError:
3168 self.del_component(self.compress_fact_phase)
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/expression.py:375, in Expression.construct(self, data)
372 try:
373 # We do not (currently) accept data for constructing Constraints
374 assert data is None
--> 375 self._construct_from_rule_using_setitem()
376 finally:
377 timer.report()
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/core/base/indexed_component.py:784, in IndexedComponent._construct_from_rule_using_setitem(self)
782 else:
783 for index in self.index_set():
--> 784 self._setitem_when_not_present(index, rule(block, index))
785 except:
786 err = sys.exc_info()[1]
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/models/properties/modular_properties/base/generic_property.py:3162, in GenericStateBlockData._compress_fact_phase.<locals>.rule_Z_phase(b, p)
3160 def rule_Z_phase(b, p):
3161 p_config = b.params.get_phase(p).config
-> 3162 return p_config.equation_of_state.compress_fact_phase(b, p)
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/models/properties/modular_properties/eos/ceos.py:551, in Cubic.compress_fact_phase(b, p)
548 A = getattr(b, cname + "_A")
549 B = getattr(b, cname + "_B")
--> 551 expr_write = CubicThermoExpressions(b)
552 if pobj.is_vapor_phase():
553 return expr_write.z_vap(eos=pobj._cubic_type, A=A[p], B=B[p])
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/models/properties/modular_properties/eos/ceos_common.py:106, in CubicThermoExpressions.__init__(self, blk)
104 def __init__(self, blk):
105 if not cubic_roots_available():
--> 106 raise RuntimeError("Cubic root external functions are not available.")
107 self.blk = blk
RuntimeError: Cubic root external functions are not available.
m.fs.R101 = GibbsReactor(
property_package=m.fs.thermo_params,
has_heat_transfer=True,
has_pressure_change=False,
)
Connecting Unit Models Using Arcs#
We have now added all the unit models we need to the flowsheet. Let us connect the unit models by defining and building each Arc:
m.fs.s01 = Arc(source=m.fs.CH4.outlet, destination=m.fs.M101.methane_feed)
m.fs.s02 = Arc(source=m.fs.H2O.outlet, destination=m.fs.M101.steam_feed)
m.fs.s03 = Arc(source=m.fs.M101.outlet, destination=m.fs.C101.inlet)
m.fs.s04 = Arc(source=m.fs.C101.outlet, destination=m.fs.H101.inlet)
m.fs.s05 = Arc(source=m.fs.H101.outlet, destination=m.fs.R101.inlet)
m.fs.s06 = Arc(source=m.fs.R101.outlet, destination=m.fs.PROD.inlet)
Finally, we can use Pyomo’s TransformationFactory to write the equality constraints on each Arc:
TransformationFactory("network.expand_arcs").apply_to(m)
Adding Expressions to Compute Operating Costs#
For this flowsheet, we are interested in computing hydrogen production in millions of pounds per year, as well as the total costs due to pressurizing, cooling, and heating utilities:
Let us add Expressions to convert the product flow from mol/s to MM lb/year of hydrogen, and to calculate the cooling, heating and compression operating costs. The total operating cost will be the sum of the costs, expressed in \$/year assuming 8000 operating hours per year (~10% downtime, which is fairly common for small scale chemical plants):
m.fs.hyd_prod = Expression(
expr=pyunits.convert(
m.fs.PROD.inlet.flow_mol[0]
* m.fs.PROD.inlet.mole_frac_comp[0, "H2"]
* m.fs.thermo_params.H2.mw, # MW defined in properties as kg/mol
to_units=pyunits.Mlb / pyunits.yr,
)
) # converting kg/s to MM lb/year
m.fs.cooling_cost = Expression(
expr=0.212e-7 * (m.fs.R101.heat_duty[0])
) # the reaction is endothermic, so R101 duty is positive
m.fs.heating_cost = Expression(
expr=2.2e-7 * m.fs.H101.heat_duty[0]
) # the stream must be heated to T_rxn, so H101 duty is positive
m.fs.compression_cost = Expression(
expr=0.12e-5 * m.fs.C101.work_isentropic[0]
) # the stream must be pressurized, so the C101 work is positive
m.fs.operating_cost = Expression(
expr=(3600 * 8000 * (m.fs.heating_cost + m.fs.cooling_cost + m.fs.compression_cost))
)
Fixing Feed Conditions#
We expect each stream to have 8 degrees of freedom, the mixer to have 0 (after both streams are accounted for), the compressor to have 2 (the pressure change and efficiency), the heater to have 1 (just the duty, since the inlet is also the outlet of M101), and the reactor to have 1 (conversion). Therefore, we have 20 degrees of freedom to specify: temperature, pressure, flow and mole fractions of all five components on both streams; compressor pressure change and efficiency; outlet heater temperature; and reactor conversion.
print(degrees_of_freedom(m))
20
Based on the literature source, we will initialize our simulation with the following values:
m.fs.CH4.outlet.mole_frac_comp[0, "CH4"].fix(1)
m.fs.CH4.outlet.mole_frac_comp[0, "H2O"].fix(1e-5)
m.fs.CH4.outlet.mole_frac_comp[0, "H2"].fix(1e-5)
m.fs.CH4.outlet.mole_frac_comp[0, "CO"].fix(1e-5)
m.fs.CH4.outlet.mole_frac_comp[0, "CO2"].fix(1e-5)
m.fs.CH4.outlet.flow_mol.fix(75 * pyunits.mol / pyunits.s)
m.fs.CH4.outlet.temperature.fix(298.15 * pyunits.K)
m.fs.CH4.outlet.pressure.fix(1e5 * pyunits.Pa)
m.fs.H2O.outlet.mole_frac_comp[0, "CH4"].fix(1e-5)
m.fs.H2O.outlet.mole_frac_comp[0, "H2O"].fix(1)
m.fs.H2O.outlet.mole_frac_comp[0, "H2"].fix(1e-5)
m.fs.H2O.outlet.mole_frac_comp[0, "CO"].fix(1e-5)
m.fs.H2O.outlet.mole_frac_comp[0, "CO2"].fix(1e-5)
m.fs.H2O.outlet.flow_mol.fix(234 * pyunits.mol / pyunits.s)
m.fs.H2O.outlet.temperature.fix(373.15 * pyunits.K)
m.fs.H2O.outlet.pressure.fix(1e5 * pyunits.Pa)
Fixing Unit Model Specifications#
For the initial problem, let us fix the compressor outlet pressure to 2 bar for now, the efficiency to 0.90 (a common assumption for compressor units), and the heater outlet temperature to 500 K. We will unfix these values later to optimize the flowsheet.
m.fs.C101.outlet.pressure.fix(pyunits.convert(2 * pyunits.bar, to_units=pyunits.Pa))
m.fs.C101.efficiency_isentropic.fix(0.90)
m.fs.H101.outlet.temperature.fix(500 * pyunits.K)
The GibbsReactor unit model calculates the amount of product and reactant based on the free energy minimization; therefore, we will specify a desired conversion and let the solver determine the reactor duty and heat transfer. For convenience, we will define the reactor conversion as the amount of methane that is converted.
m.fs.R101.conversion = Var(
initialize=0.80, bounds=(0, 1), units=pyunits.dimensionless
) # fraction
m.fs.R101.conv_constraint = Constraint(
expr=m.fs.R101.conversion
* m.fs.R101.inlet.flow_mol[0]
* m.fs.R101.inlet.mole_frac_comp[0, "CH4"]
== (
m.fs.R101.inlet.flow_mol[0] * m.fs.R101.inlet.mole_frac_comp[0, "CH4"]
- m.fs.R101.outlet.flow_mol[0] * m.fs.R101.outlet.mole_frac_comp[0, "CH4"]
)
)
m.fs.R101.conversion.fix(0.80)
For initialization, we solve a square problem (degrees of freedom = 0). Let’s check the degrees of freedom below:
print(degrees_of_freedom(m))
0
Finally, we need to initialize the each unit operation and propagate the outlet results in sequence to solve the flowsheet:
# Initialize and solve each unit operation
m.fs.CH4.initialize()
propagate_state(arc=m.fs.s01)
m.fs.H2O.initialize()
propagate_state(arc=m.fs.s02)
m.fs.M101.initialize()
propagate_state(arc=m.fs.s03)
m.fs.C101.initialize()
propagate_state(arc=m.fs.s04)
m.fs.H101.initialize()
propagate_state(arc=m.fs.s05)
m.fs.R101.initialize()
propagate_state(arc=m.fs.s06)
m.fs.PROD.initialize()
# set solver
solver = get_solver()
2023-11-02 10:26:35 [INFO] idaes.init.fs.CH4.properties: Starting initialization
2023-11-02 10:26:35 [INFO] idaes.init.fs.CH4.properties: Property initialization: optimal - Optimal Solution Found.
2023-11-02 10:26:35 [INFO] idaes.init.fs.CH4.properties: Property package initialization: optimal - Optimal Solution Found.
2023-11-02 10:26:35 [INFO] idaes.init.fs.CH4: Initialization Complete.
2023-11-02 10:26:35 [INFO] idaes.init.fs.H2O.properties: Starting initialization
2023-11-02 10:26:35 [INFO] idaes.init.fs.H2O.properties: Property initialization: optimal - Optimal Solution Found.
2023-11-02 10:26:35 [INFO] idaes.init.fs.H2O.properties: Property package initialization: optimal - Optimal Solution Found.
2023-11-02 10:26:35 [INFO] idaes.init.fs.H2O: Initialization Complete.
2023-11-02 10:26:35 [INFO] idaes.init.fs.M101.methane_feed_state: Starting initialization
2023-11-02 10:26:35 [INFO] idaes.init.fs.M101.methane_feed_state: Property initialization: optimal - Optimal Solution Found.
2023-11-02 10:26:35 [INFO] idaes.init.fs.M101.steam_feed_state: Starting initialization
2023-11-02 10:26:35 [INFO] idaes.init.fs.M101.steam_feed_state: Property initialization: optimal - Optimal Solution Found.
2023-11-02 10:26:35 [INFO] idaes.init.fs.M101.mixed_state: Starting initialization
2023-11-02 10:26:35 [INFO] idaes.init.fs.M101.mixed_state: Property initialization: optimal - Optimal Solution Found.
2023-11-02 10:26:35 [INFO] idaes.init.fs.M101.mixed_state: Property package initialization: optimal - Optimal Solution Found.
2023-11-02 10:26:35 [INFO] idaes.init.fs.M101: Initialization Complete: optimal - Optimal Solution Found
2023-11-02 10:26:35 [INFO] idaes.init.fs.C101.control_volume.properties_in: Starting initialization
2023-11-02 10:26:35 [INFO] idaes.init.fs.C101.control_volume.properties_in: Property initialization: optimal - Optimal Solution Found.
2023-11-02 10:26:35 [INFO] idaes.init.fs.C101.control_volume.properties_out: Starting initialization
2023-11-02 10:26:35 [INFO] idaes.init.fs.C101.control_volume.properties_out: Property initialization: optimal - Optimal Solution Found.
2023-11-02 10:26:35 [INFO] idaes.init.fs.C101.control_volume.properties_out: Property package initialization: optimal - Optimal Solution Found.
2023-11-02 10:26:35 [INFO] idaes.init.fs.C101.properties_isentropic: Starting initialization
2023-11-02 10:26:35 [INFO] idaes.init.fs.C101.properties_isentropic: Property initialization: optimal - Optimal Solution Found.
2023-11-02 10:26:35 [INFO] idaes.init.fs.C101.properties_isentropic: Property package initialization: optimal - Optimal Solution Found.
2023-11-02 10:26:35 [INFO] idaes.init.fs.C101: Initialization Complete: optimal - Optimal Solution Found
2023-11-02 10:26:35 [INFO] idaes.init.fs.H101.control_volume.properties_in: Starting initialization
2023-11-02 10:26:35 [INFO] idaes.init.fs.H101.control_volume.properties_in: Property initialization: optimal - Optimal Solution Found.
2023-11-02 10:26:35 [INFO] idaes.init.fs.H101.control_volume.properties_out: Starting initialization
2023-11-02 10:26:35 [INFO] idaes.init.fs.H101.control_volume.properties_out: Property initialization: optimal - Optimal Solution Found.
2023-11-02 10:26:35 [INFO] idaes.init.fs.H101.control_volume: Initialization Complete
2023-11-02 10:26:35 [INFO] idaes.init.fs.H101: Initialization Complete: optimal - Optimal Solution Found
2023-11-02 10:26:35 [INFO] idaes.init.fs.R101.control_volume.properties_in: Starting initialization
2023-11-02 10:26:35 [INFO] idaes.init.fs.R101.control_volume.properties_in: Property initialization: optimal - Optimal Solution Found.
2023-11-02 10:26:35 [INFO] idaes.init.fs.R101.control_volume.properties_out: Starting initialization
2023-11-02 10:26:35 [INFO] idaes.init.fs.R101.control_volume.properties_out: Property initialization: optimal - Optimal Solution Found.
2023-11-02 10:26:35 [INFO] idaes.init.fs.R101.control_volume: Initialization Complete
2023-11-02 10:26:36 [INFO] idaes.init.fs.R101: Initialization Complete: optimal - Optimal Solution Found
2023-11-02 10:26:36 [INFO] idaes.init.fs.PROD.properties: Starting initialization
2023-11-02 10:26:36 [INFO] idaes.init.fs.PROD.properties: Property initialization: optimal - Optimal Solution Found.
2023-11-02 10:26:36 [INFO] idaes.init.fs.PROD.properties: Property package initialization: optimal - Optimal Solution Found.
2023-11-02 10:26:36 [INFO] idaes.init.fs.PROD: Initialization Complete.
# Solve the model
results = solver.solve(m, tee=True)
Ipopt 3.13.2: nlp_scaling_method=gradient-based
tol=1e-06
max_iter=200
******************************************************************************
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...: 591
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 490
Total number of variables............................: 203
variables with only lower bounds: 13
variables with lower and upper bounds: 179
variables with only upper bounds: 0
Total number of equality constraints.................: 203
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.49e+06 1.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 0.0000000e+00 1.35e+04 2.00e-01 -1.0 3.59e+00 - 9.90e-01 9.91e-01h 1
2 0.0000000e+00 3.59e-04 9.99e+00 -1.0 3.56e+00 - 9.90e-01 1.00e+00h 1
3 0.0000000e+00 1.75e-08 8.98e+01 -1.0 2.91e-04 - 9.90e-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....: 1.1641532182693481e-10 1.7462298274040222e-08
Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00
Overall NLP error.......: 1.1641532182693481e-10 1.7462298274040222e-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.002
Total CPU secs in NLP function evaluations = 0.001
EXIT: Optimal Solution Found.
Analyze the Results of the Square Problem#
What is the total operating cost?
print(f"operating cost = ${value(m.fs.operating_cost)/1e6:0.3f} million per year")
operating cost = $39.958 million per year
For this operating cost, what conversion did we achieve of methane to hydrogen?
m.fs.R101.report()
print()
print(f"Conversion achieved = {value(m.fs.R101.conversion):.1%}")
====================================================================================
Unit : fs.R101 Time: 0.0
------------------------------------------------------------------------------------
Unit Performance
Variables:
Key : Value : Units : Fixed : Bounds
Heat Duty : 1.7819e+07 : watt : False : (None, None)
------------------------------------------------------------------------------------
Stream Table
Units Inlet Outlet
Total Molar Flowrate mole / second 309.01 429.02
Total Mole Fraction CH4 dimensionless 0.24272 0.034965
Total Mole Fraction H2O dimensionless 0.75725 0.32532
Total Mole Fraction H2 dimensionless 9.9996e-06 0.49984
Total Mole Fraction CO dimensionless 9.9996e-06 0.059609
Total Mole Fraction CO2 dimensionless 9.9996e-06 0.080265
Temperature kelvin 500.00 920.80
Pressure pascal 2.0000e+05 2.0000e+05
====================================================================================
Conversion achieved = 80.0%
Optimizing Hydrogen Production#
Now that the flowsheet has been squared and solved, we can run a small optimization problem to determine optimal conditions for producing hydrogen. Suppose we wish to find ideal conditions for the competing reactions. The GibbsReactor does not drive equilibrium forward based on temperature, so we will see small amounts of intermediate components present in the product stream. We will allow for variable reactor temperature and pressure by freeing our heater and compressor specifications, and minimize cost to achieve 90% methane conversion. Since we assume an isentopic compressor, allowing compression will heat our feed stream and reduce or eliminate the required heater duty.
Let us declare our objective function for this problem.
m.fs.objective = Objective(expr=m.fs.operating_cost)
Now, we need to add the design constraints and unfix the decision variables as we had solved a square problem until now, as well as set bounds for the design variables (reactor outlet temperature is set by state variable bounds in property package):
m.fs.R101.conversion.fix(0.90)
m.fs.C101.outlet.pressure.unfix()
m.fs.C101.outlet.pressure[0].setlb(
pyunits.convert(1 * pyunits.bar, to_units=pyunits.Pa)
) # equals inlet pressure
m.fs.C101.outlet.pressure[0].setlb(
pyunits.convert(10 * pyunits.bar, to_units=pyunits.Pa)
) # at most, pressurize to 1 bar
m.fs.H101.outlet.temperature.unfix()
m.fs.H101.heat_duty[0].setlb(
0 * pyunits.J / pyunits.s
) # ensures outlet is equal to or greater than inlet temperature
m.fs.H101.outlet.temperature[0].setub(1000 * pyunits.K) # at most, heat to 1000 K
We have now defined the optimization problem and we are now ready to solve this problem.
results = solver.solve(m, tee=True)
Ipopt 3.13.2: nlp_scaling_method=gradient-based
tol=1e-06
max_iter=200
******************************************************************************
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...: 598
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 506
Total number of variables............................: 205
variables with only lower bounds: 14
variables with lower and upper bounds: 181
variables with only upper bounds: 0
Total number of equality constraints.................: 203
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 3.9958388e+07 1.49e+06 3.46e+01 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 3.8920063e+07 1.48e+06 1.52e+03 -1.0 7.19e+06 - 3.91e-01 6.43e-03f 1
2 7.0948609e+07 1.15e+06 1.86e+06 -1.0 4.83e+06 - 1.51e-01 2.26e-01h 1
3 1.0553921e+08 5.23e+05 1.04e+07 -1.0 2.42e+06 - 3.41e-01 5.67e-01h 1
4 1.0874890e+08 1.58e+05 7.64e+06 -1.0 8.45e+05 - 7.09e-01 7.11e-01h 1
5 1.0751027e+08 1.51e+04 1.67e+06 -1.0 2.97e+05 - 9.49e-01 9.09e-01f 1
6 1.0721898e+08 5.95e+00 9.98e+03 -1.0 3.47e+04 - 9.90e-01 1.00e+00f 1
7 1.0721794e+08 3.43e-05 8.84e+01 -1.0 1.59e+02 - 9.90e-01 1.00e+00f 1
8 1.0721794e+08 8.85e-09 7.14e-01 -1.0 1.43e-02 - 9.92e-01 1.00e+00h 1
9 1.0721794e+08 1.50e-08 1.92e-06 -2.5 1.72e-02 - 1.00e+00 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 1.0721794e+08 1.02e-08 2.70e-06 -3.8 4.73e-04 - 1.00e+00 1.00e+00f 1
11 1.0721794e+08 1.02e-08 2.31e-06 -5.7 2.63e-05 - 1.00e+00 1.00e+00f 1
12 1.0721794e+08 1.49e-08 1.00e-06 -7.0 3.06e-07 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 12
(scaled) (unscaled)
Objective...............: 1.0721793780338201e+08 1.0721793780338201e+08
Dual infeasibility......: 1.0011035875619890e-06 1.0011035875619890e-06
Constraint violation....: 1.5205920451933321e-12 1.4901161193847656e-08
Complementarity.........: 9.0909090914354020e-08 9.0909090914354020e-08
Overall NLP error.......: 9.0909090914354020e-08 1.0011035875619890e-06
Number of objective function evaluations = 13
Number of objective gradient evaluations = 13
Number of equality constraint evaluations = 13
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 13
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 12
Total CPU secs in IPOPT (w/o function evaluations) = 0.004
Total CPU secs in NLP function evaluations = 0.004
EXIT: Optimal Solution Found.
print(f"operating cost = ${value(m.fs.operating_cost)/1e6:0.3f} million per year")
print()
print("Compressor results")
m.fs.C101.report()
print()
print("Heater results")
m.fs.H101.report()
print()
print("Gibbs reactor results")
m.fs.R101.report()
operating cost = $107.218 million per year
Compressor results
====================================================================================
Unit : fs.C101 Time: 0.0
------------------------------------------------------------------------------------
Unit Performance
Variables:
Key : Value : Units : Fixed : Bounds
Isentropic Efficiency : 0.90000 : dimensionless : True : (None, None)
Mechanical Work : 3.0334e+06 : watt : False : (None, None)
Pressure Change : 9.0000e+05 : pascal : False : (None, None)
Pressure Ratio : 10.000 : dimensionless : False : (None, None)
------------------------------------------------------------------------------------
Stream Table
Units Inlet Outlet
Total Molar Flowrate mole / second 309.01 309.01
Total Mole Fraction CH4 dimensionless 0.24272 0.24272
Total Mole Fraction H2O dimensionless 0.75725 0.75725
Total Mole Fraction H2 dimensionless 9.9996e-06 9.9996e-06
Total Mole Fraction CO dimensionless 9.9996e-06 9.9996e-06
Total Mole Fraction CO2 dimensionless 9.9996e-06 9.9996e-06
Temperature kelvin 353.80 619.25
Pressure pascal 1.0000e+05 1.0000e+06
====================================================================================
Heater results
====================================================================================
Unit : fs.H101 Time: 0.0
------------------------------------------------------------------------------------
Unit Performance
Variables:
Key : Value : Units : Fixed : Bounds
Heat Duty : 5.8781e-09 : watt : False : (0.0, None)
------------------------------------------------------------------------------------
Stream Table
Units Inlet Outlet
Total Molar Flowrate mole / second 309.01 309.01
Total Mole Fraction CH4 dimensionless 0.24272 0.24272
Total Mole Fraction H2O dimensionless 0.75725 0.75725
Total Mole Fraction H2 dimensionless 9.9996e-06 9.9996e-06
Total Mole Fraction CO dimensionless 9.9996e-06 9.9996e-06
Total Mole Fraction CO2 dimensionless 9.9996e-06 9.9996e-06
Temperature kelvin 619.25 619.25
Pressure pascal 1.0000e+06 1.0000e+06
====================================================================================
Gibbs reactor results
====================================================================================
Unit : fs.R101 Time: 0.0
------------------------------------------------------------------------------------
Unit Performance
Variables:
Key : Value : Units : Fixed : Bounds
Heat Duty : 2.1076e+07 : watt : False : (None, None)
------------------------------------------------------------------------------------
Stream Table
Units Inlet Outlet
Total Molar Flowrate mole / second 309.01 444.02
Total Mole Fraction CH4 dimensionless 0.24272 0.016892
Total Mole Fraction H2O dimensionless 0.75725 0.31609
Total Mole Fraction H2 dimensionless 9.9996e-06 0.51498
Total Mole Fraction CO dimensionless 9.9996e-06 0.093140
Total Mole Fraction CO2 dimensionless 9.9996e-06 0.058900
Temperature kelvin 619.25 1087.4
Pressure pascal 1.0000e+06 1.0000e+06
====================================================================================
Display optimal values for the decision variables and design variables:
print("Optimal Values")
print()
print(f"C101 outlet pressure = {value(m.fs.C101.outlet.pressure[0])/1E6:0.3f} MPa")
print()
print(f"C101 outlet temperature = {value(m.fs.C101.outlet.temperature[0]):0.3f} K")
print()
print(f"H101 outlet temperature = {value(m.fs.H101.outlet.temperature[0]):0.3f} K")
print()
print(f"R101 outlet temperature = {value(m.fs.R101.outlet.temperature[0]):0.3f} K")
print()
print(f"Hydrogen produced = {value(m.fs.hyd_prod):0.3f} MM lb/year")
print()
print(f"Conversion achieved = {value(m.fs.R101.conversion):.1%}")
Optimal Values
C101 outlet pressure = 1.000 MPa
C101 outlet temperature = 619.248 K
H101 outlet temperature = 619.248 K
R101 outlet temperature = 1087.385 K
Hydrogen produced = 32.070 MM lb/year
Conversion achieved = 90.0%