IDAES / idaes-pse

The IDAES Process Systems Engineering Framework
https://idaes-pse.readthedocs.io/
Other
216 stars 232 forks source link

Handling inert species in Gibbs Reactor #1392

Closed JavalVyas2000 closed 5 months ago

JavalVyas2000 commented 5 months ago

The inert_species_balance function in gibbs reactor does not create the balance for all the inert species, but only for the ones that are linearly independent with element balance as per the logic in the function. This increases the degrees of freedom as the variables are created for all the species including inerts but there are no constraints addressing the linearly independent with element balance.

The previous version of gibbs reactor use to create the constraints for all the inert species commit id - 119656fcfc9d249d844f680b14c74da991e6f9dd

To give an example of this, look at the below example

from idaes.models.unit_models import GibbsReactor
import pyomo.environ as pyo
from idaes.models.properties.modular_properties.base.generic_property import (
    GenericParameterBlock,
)
from idaes.models.properties.modular_properties.base.generic_reaction import (
    GenericReactionParameterBlock,
)
from idaes.models_extra.power_generation.properties.natural_gas_PR import (
    get_prop,
    get_rxn,
)

m=pyo.ConcreteModel()
m.fs=FlowsheetBlock(dynamic=False)
NG_config = get_prop(
        components=[
            "H2",
            "CO",
            "H2O",
            "CO2",
            "CH4",
            "C2H6",
            "C3H8",
            "C4H10",
            "N2",
            "O2",
            "Ar",
        ]
    )
m.fs.NG_props = GenericParameterBlock(**NG_config)
m.fs.gr=GibbsReactor(has_heat_transfer=True,
        has_pressure_change=True,
        inert_species=["N2", "Ar", "O2"],
        property_package=m.fs.NG_props)

### This is the inert_species_balance function
inert_species=["N2", "Ar", "O2"]
for j in inert_species:
    cv = m.fs.gr.control_volume
    e_comp = cv.properties_out[0].config.parameters.element_comp

    # Check for linear dependence with element balances
    # If an inert species is the only source of element e,
    # the inert species balance would be linearly dependent on the
    # element balance for e.
    dependent = True

    if len(m.fs.gr.control_volume.properties_in.phase_list) > 1:
        # Multiple phases avoid linear dependency
        dependent = False
    else:
        for e in m.fs.gr.config.property_package.element_list:
            if e_comp[j][e] == 0:
                # Element e not in component j, no effect
                continue
            else:
                for i in m.fs.gr.control_volume.properties_in.component_list:
                    if i == j:
                        continue
                    else:
                        # If comp j shares element e with comp i
                        # cannot be linearly dependent
                        if e_comp[i][e] != 0:
                            print(j,e)
                            dependent = False
    print(j,dependent)
    if (
        not dependent
        and ('Vap', j) in m.fs.gr.control_volume.properties_in.phase_component_set
    ):
        print(j)
        # return 0 == (
        #     cv.properties_in[t].get_material_flow_terms(p, j)
        #     - cv.properties_out[t].get_material_flow_terms(p, j)
        # )
    else:
        print(j)
        # return Constraint.Skip

This gives shows that the model does not create constraints for N2 and Ar as they are linearly independent, but creates constraints for O2 as it has O. This makes the dof to increase by 2 of the flowsheet for each Gibbs Reactor.

andrewlee94 commented 5 months ago

So, my first reaction here is that this is user error. How can O2 be inert if you also have CO2, CO, H2O, and hydrocarbons? What happens if O2 is removed from the inert species list?

JavalVyas2000 commented 5 months ago

I checked the DOF after removing O2 and it was the same.

JavalVyas2000 commented 5 months ago

Another possible issue I found with this is that the Gibbs Reactor makes the lagrange_mult variable for all the elements in the property package but does not consider if the element is linearly independent. As the linearly independent elements have the lagrange_mult as a free variable. This made the DOF to be 0 in this case.