Gibbs Reactor Modeling for CLHP Process – Multiphase Property Package

[quddus19]

Hello,

I’m currently working on modeling a Chemical Looping Hydrogen Production (CLHP) process using IDAES. While I am aware of the counter-current moving bed example, my current implementation is more aligned with the Aspen Plus configuration, where the system consists of a series of five Gibbs reactors.
To explore better optimization capabilities, I’m trying to replicate this setup in IDAES using Gibbs reactors. However, I'm encountering some confusion, especially related to property packages. Since the Gibbs reactor in IDAES requires a single property package, I attempted to include both gas and solid phases in one package. This led to the following error:

KeyError: "Index '('Vap', 'Fe')' is not valid for indexed component 'fs.reactor.control_volume.properties_in[0.0].flow_mol_phase_comp'"

I haven't yet found a working example where a single property package handles multiphase (Vap + Sol) components for a Gibbs reactor.
Are there any existing examples that demonstrate this, or would you have any recommendations on how to properly set up a unified multiphase property package?
Also, I understand that in general, solids do not require an equation of state (EOS). Is it still mandatory to provide one in IDAES for the solid phase, or is there a way around it?

Thank you for your help!

[dallan-keylogic]

Can you please provide the entire stack trace, not just the final error message? While theoretically the Gibbs reactor should work for multiphase systems, it has not been tested.

IDAES (or, more precisely, the modular property framework) uses the term "equation of state" more broadly than is traditional in the thermodynamic literature. For example, the Ideal equation of state treats gas phases as ideal gases but treats liquid and solid phases as ideal mixtures. Activity coefficient models are also considered "equations of state", with a partial implementation of eNRTL existing in the IDAES repo.

[quddus19]

Thank you again for your response, and I really appreciate the explanation about the equation of state.
Below is my entire stack trace.

ERROR: Rule failed for Expression
'fs.reactor.control_volume.elemental_flow_in' with index (0.0, 'Vap', 'C'):
KeyError: "Index '('Vap', 'Fe')' is not valid for indexed component
'fs.reactor.control_volume.properties_in[0.0].flow_mol_phase_comp'"
ERROR: Constructing component 'fs.reactor.control_volume.elemental_flow_in'
from data=None failed:
        KeyError: "Index '('Vap', 'Fe')' is not valid for indexed component
        'fs.reactor.control_volume.properties_in[0.0].flow_mol_phase_comp'"
2025-04-24 19:44:27 [ERROR] idaes.core.base.process_block: Failure in build: fs.reactor
Traceback (most recent call last):
  File "C:\Users\Bang\anaconda3\Lib\site-packages\idaes\core\base\process_block.py", line 41, in _rule_default
    b.build()
  File "C:\Users\Bang\AppData\Local\Temp\ipykernel_88952\429995452.py", line 363, in build
    self.control_volume.add_total_element_balances()
  File "C:\Users\Bang\anaconda3\Lib\site-packages\idaes\core\base\control_volume0d.py", line 1037, in add_total_element_balances
    @self.Expression(
     ^^^^^^^^^^^^^^^^
  File "C:\Users\Bang\anaconda3\Lib\site-packages\pyomo\core\base\block.py", line 77, in __call__
    setattr(
  File "C:\Users\Bang\anaconda3\Lib\site-packages\pyomo\core\base\block.py", line 571, in __setattr__
    self.add_component(name, val)
  File "C:\Users\Bang\anaconda3\Lib\site-packages\pyomo\core\base\block.py", line 1101, in add_component
    val.construct(data)
  File "C:\Users\Bang\anaconda3\Lib\site-packages\pyomo\core\base\expression.py", line 383, in construct
    self._construct_from_rule_using_setitem()
  File "C:\Users\Bang\anaconda3\Lib\site-packages\pyomo\core\base\indexed_component.py", line 788, in _construct_from_rule_using_setitem
    self._setitem_when_not_present(index, rule(block, index))
                                          ^^^^^^^^^^^^^^^^^^
  File "C:\Users\Bang\anaconda3\Lib\site-packages\pyomo\core\base\initializer.py", line 347, in __call__
    return self._fcn(parent, *idx)
           ^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\Bang\anaconda3\Lib\site-packages\idaes\core\base\control_volume0d.py", line 1044, in elemental_flow_in
    return sum(
           ^^^^
  File "C:\Users\Bang\anaconda3\Lib\site-packages\idaes\core\base\control_volume0d.py", line 1046, in <genexpr>
    * b.properties_in[t].get_material_flow_terms(p, j)
      ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\Bang\anaconda3\Lib\site-packages\idaes\models\properties\modular_properties\state_definitions\FTPx.py", line 270, in get_material_flow_terms_FTPx
    return b.flow_mol_phase_comp[p, j]
           ~~~~~~~~~~~~~~~~~~~~~^^^^^^
  File "C:\Users\Bang\anaconda3\Lib\site-packages\pyomo\core\base\indexed_component.py", line 648, in __getitem__
    validated_index = self._validate_index(index)
                      ^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\Bang\anaconda3\Lib\site-packages\pyomo\core\base\indexed_component.py", line 870, in _validate_index
    raise KeyError(
KeyError: "Index '('Vap', 'Fe')' is not valid for indexed component 'fs.reactor.control_volume.properties_in[0.0].flow_mol_phase_comp'"
ERROR: Constructing component 'fs.reactor' from data=None failed:
        KeyError: "Index '('Vap', 'Fe')' is not valid for indexed component
        'fs.reactor.control_volume.properties_in[0.0].flow_mol_phase_comp'"
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
Cell In[8], line 18
     16 # 3. 유닛 구성
     17 m.fs.feed = Feed(property_package=m.fs.thermo_props)
---> 18 m.fs.reactor = GibbsReactor2(
     19     property_package=m.fs.thermo_props,
     20     has_heat_transfer=False,
     21     has_pressure_change=False,
     22 )
     23 m.fs.product = Product(property_package=m.fs.thermo_props)
     25 # 4. 연결

File ~\anaconda3\Lib\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 ~\anaconda3\Lib\site-packages\pyomo\core\base\block.py:1101, in BlockData.add_component(self, name, val)
   1093     logger.debug(
   1094         "Constructing %s '%s' on %s from data=%s",
   1095         val.__class__.__name__,
   (...)
   1098         str(data),
   1099     )
   1100 try:
-> 1101     val.construct(data)
   1102 except:
   1103     err = sys.exc_info()[1]

File ~\anaconda3\Lib\site-packages\pyomo\core\base\block.py:2158, in Block.construct(self, data)
   2156                     obj.construct(data.get(name, None))
   2157         # Trigger the (normal) initialization of the block
-> 2158         self._getitem_when_not_present(_idx)
   2159 finally:
   2160     # We must allow that id(self) may no longer be in
   2161     # _BlockConstruction.data, as _getitem_when_not_present will
   2162     # have already removed the entry for scalar blocks (as the
   2163     # BlockData and the Block component are the same object)
   2164     if data is not None:

File ~\anaconda3\Lib\site-packages\pyomo\core\base\block.py:2073, in Block._getitem_when_not_present(self, idx)
   2070     data = None
   2072 try:
-> 2073     obj = self._rule(_block, idx)
   2074     # If the user returns a block, transfer over everything
   2075     # they defined into the empty one we created.  We do
   2076     # this inside the try block so that any abstract
   2077     # components declared by the rule have the opportunity
   2078     # to be initialized with data from
   2079     # _BlockConstruction.data as they are transferred over.
   2080     if obj is not _block and isinstance(obj, BlockData):

File ~\anaconda3\Lib\site-packages\pyomo\core\base\initializer.py:349, in IndexedCallInitializer.__call__(self, parent, idx)
    347     return self._fcn(parent, *idx)
    348 else:
--> 349     return self._fcn(parent, idx)

File ~\anaconda3\Lib\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}")

Cell In[1], line 363, in GibbsReactorData.build(self)
    357     if i not in self.control_volume.properties_in.component_list:
    358         raise ConfigurationError(
    359             "{} invalid component in inert_species argument. {} is "
    360             "not in the property package component list.".format(self.name, i)
    361         )
--> 363 self.control_volume.add_total_element_balances()
    365 self.control_volume.add_energy_balances(
    366     balance_type=self.config.energy_balance_type,
    367     has_heat_transfer=self.config.has_heat_transfer,
    368 )
    370 self.control_volume.add_momentum_balances(
    371     balance_type=self.config.momentum_balance_type,
    372     has_pressure_change=self.config.has_pressure_change,
    373 )

File ~\anaconda3\Lib\site-packages\idaes\core\base\control_volume0d.py:1037, in ControlVolume0DBlockData.add_total_element_balances(self, has_rate_reactions, has_equilibrium_reactions, has_phase_equilibrium, has_mass_transfer, custom_elemental_term)
   1031     else:
   1032         raise BalanceTypeNotSupportedError(
   1033             "{} property package MaterialFlowBasis == 'other'. Cannot "
   1034             "automatically generate elemental balances.".format(self.name)
   1035         )
-> 1037 @self.Expression(
   1038     self.flowsheet().time,
   1039     phase_list,
   1040     self.config.property_package.element_list,
   1041     doc="Inlet elemental flow terms",
   1042 )
   1043 def elemental_flow_in(b, t, p, e):
   1044     return sum(
   1045         conv_factor(b, t, j)
   1046         * b.properties_in[t].get_material_flow_terms(p, j)
   1047         * b.properties_out[t].params.element_comp[j][e]
   1048         for j in component_list
   1049     )
   1051 @self.Expression(
   1052     self.flowsheet().time,
   1053     phase_list,
   (...)
   1056 )
   1057 def elemental_flow_out(b, t, p, e):

File ~\anaconda3\Lib\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 ~\anaconda3\Lib\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 ~\anaconda3\Lib\site-packages\pyomo\core\base\block.py:1101, in BlockData.add_component(self, name, val)
   1093     logger.debug(
   1094         "Constructing %s '%s' on %s from data=%s",
   1095         val.__class__.__name__,
   (...)
   1098         str(data),
   1099     )
   1100 try:
-> 1101     val.construct(data)
   1102 except:
   1103     err = sys.exc_info()[1]

File ~\anaconda3\Lib\site-packages\pyomo\core\base\expression.py:383, in Expression.construct(self, data)
    380 try:
    381     # We do not (currently) accept data for constructing Constraints
    382     assert data is None
--> 383     self._construct_from_rule_using_setitem()
    384 finally:
    385     timer.report()

File ~\anaconda3\Lib\site-packages\pyomo\core\base\indexed_component.py:788, in IndexedComponent._construct_from_rule_using_setitem(self)
    786     else:
    787         for index in self.index_set():
--> 788             self._setitem_when_not_present(index, rule(block, index))
    789 except:
    790     err = sys.exc_info()[1]

File ~\anaconda3\Lib\site-packages\pyomo\core\base\initializer.py:347, in IndexedCallInitializer.__call__(self, parent, idx)
    341 def __call__(self, parent, idx):
    342     # Note: this is called by a component using data from a Set (so
    343     # any tuple-like type should have already been checked and
    344     # converted to a tuple; or flattening is turned off and it is
    345     # the user's responsibility to sort things out.
    346     if idx.__class__ is tuple:
--> 347         return self._fcn(parent, *idx)
    348     else:
    349         return self._fcn(parent, idx)

File ~\anaconda3\Lib\site-packages\idaes\core\base\control_volume0d.py:1044, in ControlVolume0DBlockData.add_total_element_balances.<locals>.elemental_flow_in(b, t, p, e)
   1037 @self.Expression(
   1038     self.flowsheet().time,
   1039     phase_list,
   (...)
   1042 )
   1043 def elemental_flow_in(b, t, p, e):
-> 1044     return sum(
   1045         conv_factor(b, t, j)
   1046         * b.properties_in[t].get_material_flow_terms(p, j)
   1047         * b.properties_out[t].params.element_comp[j][e]
   1048         for j in component_list
   1049     )

File ~\anaconda3\Lib\site-packages\idaes\core\base\control_volume0d.py:1046, in <genexpr>(.0)
   1037 @self.Expression(
   1038     self.flowsheet().time,
   1039     phase_list,
   (...)
   1042 )
   1043 def elemental_flow_in(b, t, p, e):
   1044     return sum(
   1045         conv_factor(b, t, j)
-> 1046         * b.properties_in[t].get_material_flow_terms(p, j)
   1047         * b.properties_out[t].params.element_comp[j][e]
   1048         for j in component_list
   1049     )

File ~\anaconda3\Lib\site-packages\idaes\models\properties\modular_properties\state_definitions\FTPx.py:270, in define_state.<locals>.get_material_flow_terms_FTPx(b, p, j)
    268 def get_material_flow_terms_FTPx(b, p, j):
    269     """Create material flow terms for control volume."""
--> 270     return b.flow_mol_phase_comp[p, j]

File ~\anaconda3\Lib\site-packages\pyomo\core\base\indexed_component.py:648, in IndexedComponent.__getitem__(self, index)
    646 if isinstance(index, EXPR.GetItemExpression):
    647     return index
--> 648 validated_index = self._validate_index(index)
    649 if validated_index is not index:
    650     index = validated_index

File ~\anaconda3\Lib\site-packages\pyomo\core\base\indexed_component.py:870, in IndexedComponent._validate_index(self, idx)
    863     raise KeyError(
    864         "Cannot treat the scalar component '%s' "
    865         "as an indexed component" % (self.name,)
    866     )
    867 #
    868 # Raise an exception
    869 #
--> 870 raise KeyError(
    871     "Index '%s' is not valid for indexed component '%s'"
    872     % (normalized_idx, self.name)
    873 )

KeyError: "Index '('Vap', 'Fe')' is not valid for indexed component 'fs.reactor.control_volume.properties_in[0.0].flow_mol_phase_comp'"  

Would it be alright to share my property package here? It's not perfect, but I hope it’s useful.

[dallan-keylogic]

Yes, feel free to share your property package. I'll take a look at the stack trace.

[quddus19]

Here's my property package code. I adapted much of it from the natural gas property example.

# Import Python libraries
import logging
import copy
import enum

# Import Pyomo units
import pyomo.environ as pyo
from pyomo.environ import units as pyunits

# Import IDAES cores
from idaes.core import VaporPhase, SolidPhase, Component, PhaseType    
from idaes.models.properties.modular_properties.state_definitions import FTPx
from idaes.models.properties.modular_properties.eos.ceos import Cubic, CubicType
from idaes.models.properties.modular_properties.eos.ideal import Ideal
from idaes.models.properties.modular_properties.phase_equil.forms import (
    log_fugacity,
)
from idaes.models.properties.modular_properties.phase_equil import SmoothVLE
from idaes.models.properties.modular_properties.pure import (
    NIST,
    RPP4,
    RPP5,
    ChapmanEnskogLennardJones,
    Eucken,
)
from idaes.models.properties.modular_properties.base.generic_reaction import (
    ConcentrationForm,
)
from idaes.models.properties.modular_properties.transport_properties import (
    ViscosityWilke,
    ThermalConductivityWMS,
    NoMethod,
)
from idaes.models.properties.modular_properties.transport_properties.viscosity_wilke import (
    wilke_phi_ij_callback,
)

from idaes.models.properties.modular_properties.reactions.dh_rxn import constant_dh_rxn
from idaes.models.properties.modular_properties.reactions.rate_constant import arrhenius
from idaes.models.properties.modular_properties.reactions.rate_forms import (
    power_law_rate,
)
from idaes.core.util.exceptions import ConfigurationError

# Set up logger
_log = logging.getLogger(__name__)


class EosType(enum.Enum):
    PR = 1
    IDEAL = 2


# Property Sources

# Source: NIST webbook
# Properties: Heat capacity coefficients for all species except ethane,
# propane, and butane. Reference enthalpies and entropies for all species.

# Source: The Properties of Gases and Liquids (1987)
# 4th edition, Chemical Engineering Series - Robert C. Reid
# Properties: Critical temperatures and pressures. Omega.
# Heat capacity coefficients for ethane, propane, and butane.


_phase_dicts_ideal = {
    "Vap": {
        "type": VaporPhase,
        "equation_of_state": Ideal,
        "visc_d_phase": ViscosityWilke,
        "transport_property_options": {
            "viscosity_phi_ij_callback": wilke_phi_ij_callback,
        },
        "therm_cond_phase": ThermalConductivityWMS,
    },
    "Sol": {
        "type": SolidPhase,
        #"equation_of_state": Ideal,
    },
}

_component_params = {
    "H2": {
        "type": Component,
        "valid_phase_types": [PhaseType.vaporPhase],
        "elemental_composition": {"H": 2},
        "enth_mol_ig_comp": NIST,
        "entr_mol_ig_comp": NIST,            
        "cp_mol_ig_comp": NIST,
        "visc_d_phase_comp": {"Vap": ChapmanEnskogLennardJones},
        "therm_cond_phase_comp": {"Vap": Eucken},
        "parameter_data": {
            "mw": (0.0020159, pyunits.kg / pyunits.mol),
            "pressure_crit": (13e5, pyunits.Pa),
            "temperature_crit": (33.2, pyunits.K),
            "omega": -0.218,
            "cp_mol_ig_comp_coeff": {
                "A": 33.066178,
                "B": -11.363417,
                "C": 11.432816,
                "D": -2.772874,
                "E": -0.158558,
                "F": -9.980797,
                "G": 172.707974,
                "H": 0.0,
            },
            "lennard_jones_sigma": (2.826, pyunits.angstrom),
            "lennard_jones_epsilon_reduced": (59.7, pyunits.K),
            "f_int_eucken": 1,
        },
    },
    "CO": {
        "type": Component,
        "valid_phase_types": [PhaseType.vaporPhase],
        "elemental_composition": {"C": 1, "O": 1},
        "enth_mol_ig_comp": NIST,
        "entr_mol_ig_comp": NIST,
        "cp_mol_ig_comp": NIST,
        "visc_d_phase_comp": {"Vap": ChapmanEnskogLennardJones},
        "therm_cond_phase_comp": {"Vap": Eucken},
        "parameter_data": {
            "mw": (0.0280101, pyunits.kg / pyunits.mol),
            "pressure_crit": (35e5, pyunits.Pa),
            "temperature_crit": (132.9, pyunits.K),
            "omega": 0.066,
            "cp_mol_ig_comp_coeff": {
                "A": 25.56759,
                "B": 6.09613,
                "C": 4.054656,
                "D": -2.671301,
                "E": 0.131021,
                "F": -118.0089,
                "G": 227.3665,
                "H": -110.5271,
            },
            "lennard_jones_sigma": (3.690, pyunits.angstrom),
            "lennard_jones_epsilon_reduced": (91.7, pyunits.K),
            "f_int_eucken": 1,
        },
    },
    "H2O": {
        "type": Component,
        "valid_phase_types": [PhaseType.vaporPhase],
        "elemental_composition": {"H": 2, "O": 1},
        "enth_mol_ig_comp": NIST,
        "entr_mol_ig_comp": NIST,
        "cp_mol_ig_comp": NIST,
        "pressure_sat_comp": NIST,
        "parameter_data": {
            "mw": (0.01801528, pyunits.kg / pyunits.mol),
            "pressure_crit": (221.2e5, pyunits.Pa),
            "temperature_crit": (647.3, pyunits.K),
            "omega": 0.344,
            "cp_mol_ig_comp_coeff": {
                "A": 30.092,
                "B": 6.832514,
                "C": 6.793435,
                "D": -2.53448,
                "E": 0.082139,
                "F": -250.881,
                "G": 223.3967,
                "H": -241.8264,
            },
            "pressure_sat_comp_coeff": {  # NIST <- Stull 1947
                "A": 4.6543,
                "B": 1435.264,
                "C": -64.848,
            },
            "lennard_jones_sigma": (2.641, pyunits.angstrom),
            "lennard_jones_epsilon_reduced": (809.1, pyunits.K),
            "f_int_eucken": 1,
        },
    },
    "CO2": {
        "type": Component,
        "valid_phase_types": [PhaseType.vaporPhase],
        "visc_d_phase_comp": {"Vap": ChapmanEnskogLennardJones},
        "therm_cond_phase_comp": {"Vap": Eucken},
        "elemental_composition": {"C": 1, "O": 2},
        "enth_mol_ig_comp": NIST,
        "entr_mol_ig_comp": NIST,
        "cp_mol_ig_comp": NIST,
        "parameter_data": {
            "mw": (0.04401, pyunits.kg / pyunits.mol),
            "pressure_crit": (73.8e5, pyunits.Pa),
            "temperature_crit": (304.1, pyunits.K),
            "omega": 0.239,
            "cp_mol_ig_comp_coeff": {
                "A": 24.99735,
                "B": 55.18696,
                "C": -33.69137,
                "D": 7.948387,
                "E": -0.136638,
                "F": -403.6075,
                "G": 228.2431,
                "H": -393.5224,
            },
            "lennard_jones_sigma": (3.941, pyunits.angstrom),
            "lennard_jones_epsilon_reduced": (195.2, pyunits.K),
            "f_int_eucken": 1,
        },
    },
    "O2": {
        "type": Component,
        "valid_phase_types": [PhaseType.vaporPhase],
        "visc_d_phase_comp": {"Vap": ChapmanEnskogLennardJones},
        "therm_cond_phase_comp": {"Vap": Eucken},
        "elemental_composition": {"O": 2},
        "enth_mol_ig_comp": NIST,
        "entr_mol_ig_comp": NIST,
        "cp_mol_ig_comp": NIST,
        "parameter_data": {
            "mw": (0.031998, pyunits.kg / pyunits.mol),
            "pressure_crit": (50.4e5, pyunits.Pa),
            "temperature_crit": (154.6, pyunits.K),
            "omega": 0.025,
            "cp_mol_ig_comp_coeff": {
                "A": 30.03235,
                "B": 8.772972,
                "C": -3.988133,
                "D": 0.788313,
                "E": -0.741599,
                "F": -11.32468,
                "G": 236.1663,
                "H": 0.0,
            },
            "lennard_jones_sigma": (3.467, pyunits.angstrom),
            "lennard_jones_epsilon_reduced": (106.7, pyunits.K),
            "f_int_eucken": 1,
        },
    },
    
    "CH4": {
        "type": Component,
        "valid_phase_types": [PhaseType.vaporPhase],
        "elemental_composition": {"C": 1, "H": 4},
        "enth_mol_ig_comp": NIST,
        "entr_mol_ig_comp": NIST,
        "cp_mol_ig_comp": NIST,
        "visc_d_phase_comp": {"Vap": ChapmanEnskogLennardJones},
        "therm_cond_phase_comp": {"Vap": Eucken},
        "parameter_data": {
            "mw": (0.0160425, pyunits.kg / pyunits.mol),
            "pressure_crit": (46e5, pyunits.Pa),
            "temperature_crit": (190.4, pyunits.K),
            "omega": 0.011,
            "cp_mol_ig_comp_coeff": {
                "A": -0.703029,
                "B": 108.4773,
                "C": -42.52157,
                "D": 5.862788,
                "E": 0.678565,
                "F": -76.84376,
                "G": 158.7163,
                "H": -74.8731,
            },
            "lennard_jones_sigma": (3.758, pyunits.angstrom),
            "lennard_jones_epsilon_reduced": (148.6, pyunits.K),
            "f_int_eucken": 1,
        },
    },
}    

_component_params_solid = {
    "Fe": {
        "type": Component,
        "valid_phase_types": [PhaseType.solidPhase],
        "elemental_composition": {"Fe": 1}, 
        "enth_mol_sol_comp": SolidNIST,
        "cp_mol_sol_comp": SolidNIST,
        "entr_mol_sol_comp": SolidNIST,
        "dens_mol_sol_comp": SolidNIST,
        "parameter_data": {
            "mw": (0.055845, pyunits.kg / pyunits.mol),
            "enth_mol_form_ref": (0, pyunits.J / pyunits.mol),
            "dens_mol": (70000, pyunits.mol / pyunits.m**3),
            "cp_mol_ig_comp_coeff": {
                "A": 23.97449, "B": 8.367750, "C": 0.000277, "D": -0.000088,
                "E": 0, "F": 0, "G": 0, "H": 0
            },
        },
    },
    "FeO": {
        "type": Component,
        "valid_phase_types": [PhaseType.solidPhase],
        "elemental_composition": {"Fe": 1, "O": 1},
        "enth_mol_sol_comp": SolidNIST,
        "cp_mol_sol_comp": SolidNIST,
        "entr_mol_sol_comp": SolidNIST,
        "dens_mol_sol_comp": SolidNIST,
        "parameter_data": {
            "mw": (0.071844, pyunits.kg / pyunits.mol),
            "enth_mol_form_ref": (-272000, pyunits.J / pyunits.mol),
            "dens_mol": (60000, pyunits.mol / pyunits.m**3),
            "cp_mol_ig_comp_coeff": {
                "A": 45.75120, "B": 18.78553, "C": -5.952201, "D": 0.852779,
                "E": 0, "F": 0, "G": 0, "H": 0
            },
        },
    },
    "Fe3O4": {
        "type": Component,
        "valid_phase_types": [PhaseType.solidPhase],
        "elemental_composition": {"Fe": 3, "O": 4}, 
        "enth_mol_sol_comp": SolidNIST,
        "cp_mol_sol_comp": SolidNIST,
        "entr_mol_sol_comp": SolidNIST,
        "dens_mol_sol_comp": SolidNIST,
        "parameter_data": {
            "mw": (0.231533, pyunits.kg / pyunits.mol),
            "enth_mol_form_ref": (-1118000, pyunits.J / pyunits.mol),
            "dens_mol": (25000, pyunits.mol / pyunits.m**3),
            "cp_mol_ig_comp_coeff": {
                "A": 104.2096, "B": 178.5108, "C": -68.2520, "D": 9.1979,
                "E": 0, "F": 0, "G": 0, "H": 0
            },
        },
    },
    "Fe2O3": {
        "type": Component,
        "valid_phase_types": [PhaseType.solidPhase],
        "elemental_composition": {"Fe": 2, "O": 3},
        "enth_mol_sol_comp": SolidNIST,
        "cp_mol_sol_comp": SolidNIST,
        "entr_mol_sol_comp": SolidNIST,
        "dens_mol_sol_comp": SolidNIST,
        "parameter_data": {
            "mw": (0.15969, pyunits.kg / pyunits.mol),
            "enth_mol_form_ref": (-824200, pyunits.J / pyunits.mol),
            "dens_mol": (30000, pyunits.mol / pyunits.m**3),
            "cp_mol_ig_comp_coeff": {
                "A": 110.9362, "B": 32.04714, "C": -9.192333, "D": 0.901506,
                "E": 0, "F": 0, "G": 0, "H": 0
            },
        },
    },
}


def get_prop(components=None, phases=["Vap", "Sol"], eos=EosType.IDEAL, scaled=False):
    if components is None:
        components = list(_component_params.keys()) + list(_component_params_solid.keys())

    configuration = {
        "components": {},
        "parameter_data": {},
        "phases": {
            "Vap": copy.deepcopy(_phase_dicts_ideal["Vap"]),
            "Sol": copy.deepcopy(_phase_dicts_ideal["Sol"]),
        },
        "base_units": {
            "time": pyunits.s,
            "length": pyunits.m,
            "mass": pyunits.kg,
            "amount": pyunits.mol,
            "temperature": pyunits.K,
        },
        "state_definition": FTPx,
        "state_bounds": {
            "flow_mol": (0, 8000, 50000, pyunits.mol / pyunits.s),
            "temperature": (273.15, 500, 2500, pyunits.K),
            "pressure": (5e4, 1.3e5, 1e8, pyunits.Pa),
        },
        "pressure_ref": (101325, pyunits.Pa),
        "temperature_ref": (298.15, pyunits.K),
        "parameter_data": {},
    }

    c = configuration["components"]
    
    for comp in components:
        if comp in _component_params:
            c[comp] = copy.deepcopy(_component_params[comp])
            if PhaseType.vaporPhase in c[comp]["valid_phase_types"]:
                configuration["phases"]["Vap"] = copy.deepcopy(_phase_dicts_ideal["Vap"])
        elif comp in _component_params_solid:
            c[comp] = copy.deepcopy(_component_params_solid[comp])
            if PhaseType.solidPhase in c[comp]["valid_phase_types"]:
                configuration["phases"]["Sol"] = copy.deepcopy(_phase_dicts_ideal["Sol"])
        else:
            raise ConfigurationError(f"Component '{comp}' not found.")

    return configuration

[dallan-keylogic]

You first need to uncomment Ideal for the solid phase.

Next, there's an error in get_prop:

    if components is None:
        components = list(_component_params.keys()) + list(_component_params_solid.keys())

This is incorrect, but it is incorrect in natural_gas_PR.py as well. It should read

    if components is None:
        components = list(comp for comp in _component_params.keys()) + list(comp for comp in _component_params_solid.keys())

Observe:

In [1]: tmp = {"A":1, "B":2}
In [2]: tmp2 = [tmp.keys()]
In [3]: tmp2
Out[3]: [dict_keys(['A', 'B'])]
In [4]: tmp3 = [key for key in tmp.keys()]
In [5]: tmp3
Out[5]: ['A', 'B']

However, the main problem isn't either of those things, but a bug in the ControlVolume0D. When the elemental_flow_in Expression is defined, it incorrectly sums over all components, instead of summing over the components in the phase.

The fix is straightforward, but testing, etc. adds additional overhead. I've opened an issue #1615 , but I make no promises about when it will be fixed. Feel free to attempt a fix yourself and open a PR.

[quddus19]

Thank you for your helpful comment. I appreciate the clarification, and I'll do my best to address the main problem regarding the incorrect component summation in elemental_flow_in. I'll look into a fix and see if I can open a PR soon.