Error Handling with Gibbs Reactor for Methanol Synthesis in IDAES

[panyaoyu]

I'm using a Gibbs reactor in the flowsheet for the synthesis of methanol, which may involve a liquid phase during the reaction process. Therefore, I need to include both vapor and liquid phases in the Property Package.

However, my reactants include components like H2, CO, and CO2, which do not have a liquid phase.

This leads to errors similar to the following:

KeyError: "Index ('Liq', 'CO') is not valid for indexed component 'fs.R101.control_volume.properties_in[0.0].flow_mol_phase_comp'"
ERROR: Constructing component 'fs.R101' from data=None failed:
        KeyError: "Index ('Liq', 'CO') is not valid for indexed component 'fs.R101.control_volume.properties_in[0.0].flow_mol_phase_comp'"

How to resolve this issue?

Thanks a lot :)

[andrewlee94]

Could you please provide the full traceback? We need to see exactly where the exception is coming from to be sure about the problem.

Additionally, it would also help if you could provide more details about your model:

1. How are you defining the property package? Is this done using the modular properties or a custom property package?
2. What approach did you use to indicate that some species do not existing in the liquid phase? It appears that you have at least partially achieved this as ("Liq", "CO2") is not a valid index, but it is possible that this was not reflected in other parts of the code.

[andrewlee94]

For the degrees of freedom issue, start by using the DiagnosticsToolbox to determine what the variables and constraints involved are.

[panyaoyu]

I think I see the cause of the first problem. Assuming you know where the core IDAES code is installed on your machine, could you try the following:

Change idaes/core/base/control_volume0d.py line 1026 to the following:

        def elemental_flow_in(b, t, p, e):
            return sum(
                conv_factor(b, t, j)
                * b.properties_in[t].get_material_flow_terms(p, j)
                * b.properties_out[t].params.element_comp[j][e]
                for j in component_list
                if (p, j) in phase_component_set  # add this line
            )

and line 1040 to:

def elemental_flow_out(b, t, p, e):
            return sum(
                conv_factor(b, t, j)
                * b.properties_out[t].get_material_flow_terms(p, j)
                * b.properties_out[t].params.element_comp[j][e]
                for j in component_list
                if (p, j) in phase_component_set  # add this line
            )

I think this will resolve this issue - it looks like the elemental balance code was not updated to allow for cases where species were not present in all phases.

Assuming this works (and you are able to test it), would you be able to open a Pull Request to make this changes (and add a test case for this too)? Unfortunately I am out on leave for the next few weeks so I will not be able to fix this for a while.

Thanks a lot! I will let you know immediately after testing :)

[panyaoyu]

I think I see the cause of the first problem. Assuming you know where the core IDAES code is installed on your machine, could you try the following:

Change idaes/core/base/control_volume0d.py line 1026 to the following:

        def elemental_flow_in(b, t, p, e):
            return sum(
                conv_factor(b, t, j)
                * b.properties_in[t].get_material_flow_terms(p, j)
                * b.properties_out[t].params.element_comp[j][e]
                for j in component_list
                if (p, j) in phase_component_set  # add this line
            )

and line 1040 to:

def elemental_flow_out(b, t, p, e):
            return sum(
                conv_factor(b, t, j)
                * b.properties_out[t].get_material_flow_terms(p, j)
                * b.properties_out[t].params.element_comp[j][e]
                for j in component_list
                if (p, j) in phase_component_set  # add this line
            )

I think this will resolve this issue - it looks like the elemental balance code was not updated to allow for cases where species were not present in all phases.

Assuming this works (and you are able to test it), would you be able to open a Pull Request to make this changes (and add a test case for this too)? Unfortunately I am out on leave for the next few weeks so I will not be able to fix this for a while.

I modified the source code of idaes according to the method you provided, and then I conducted a test. The degrees of freedom in my GibbsReactor model changed. I independently wrote a flowchart that only includes a GibbsReactor with has_heat_transfer=True, and then found that its degrees of freedom are 7, which I think is not reasonable. Because even just the FTPx of the inlet includes 5 + 1 + 1 + 1 = 8 degrees of freedom (5 components), plus the heat load should be 9 degrees of freedom. I used the display_overconstrained_set method of DiagnosticsToolbox to print the results as follows, I am very puzzled and hope to get your answer.

====================================================================================
Dulmage-Mendelsohn Over-Constrained Set

    Independent Block 0:

        Variables:

            fs.R101.control_volume.properties_out[0.0].log_mole_frac_phase_comp[Vap,CO2]
            fs.R101.control_volume.properties_out[0.0].mole_frac_phase_comp[Vap,CO2]
            fs.R101.lagrange_mult[0.0,O]
            fs.R101.lagrange_mult[0.0,C]
            fs.R101.control_volume.properties_out[0.0].mole_frac_phase_comp[Vap,CO]
            fs.R101.control_volume.properties_out[0.0].mole_frac_phase_comp[Vap,H2O]
            fs.R101.control_volume.properties_out[0.0].mole_frac_phase_comp[Vap,H2]
            fs.R101.control_volume.properties_out[0.0].mole_frac_phase_comp[Vap,CH3OH]
            fs.R101.control_volume.properties_out[0.0].flow_mol_phase[Vap]
            fs.R101.control_volume.properties_out[0.0].log_mole_frac_phase_comp[Vap,CO]
            fs.R101.lagrange_mult[0.0,H]
            fs.R101.control_volume.properties_out[0.0].log_mole_frac_phase_comp[Vap,CH3OH]
            fs.R101.control_volume.properties_out[0.0].flow_mol_phase[Liq]
            fs.R101.control_volume.properties_out[0.0].mole_frac_phase_comp[Liq,H2O]
            fs.R101.control_volume.properties_out[0.0].log_mole_frac_phase_comp[Vap,H2O]
            fs.R101.control_volume.properties_out[0.0].log_mole_frac_phase_comp[Liq,CH3OH]
            fs.R101.control_volume.properties_out[0.0]._teq[Vap,Liq]
            fs.R101.control_volume.properties_out[0.0].mole_frac_phase_comp[Liq,CH3OH]
            fs.R101.control_volume.properties_out[0.0].temperature_dew[Vap,Liq]
            fs.R101.control_volume.properties_out[0.0].log_mole_frac_phase_comp[Liq,H2O]
            fs.R101.control_volume.properties_out[0.0]._mole_frac_tdew[Vap,Liq,CO]
            fs.R101.control_volume.properties_out[0.0]._mole_frac_tdew[Vap,Liq,CO2]
            fs.R101.control_volume.properties_out[0.0].log_mole_frac_tdew[Vap,Liq,H2O]
            fs.R101.control_volume.properties_out[0.0].log_mole_frac_comp[H2O]
            fs.R101.control_volume.properties_out[0.0]._mole_frac_tdew[Vap,Liq,H2O]
            fs.R101.control_volume.properties_out[0.0]._mole_frac_tdew[Vap,Liq,H2]
            fs.R101.control_volume.properties_out[0.0]._mole_frac_tdew[Vap,Liq,CH3OH]
            fs.R101.control_volume.properties_out[0.0].log_mole_frac_tdew[Vap,Liq,CH3OH]
            fs.R101.control_volume.properties_out[0.0].log_mole_frac_comp[CH3OH]
            fs.R101.control_volume.properties_out[0.0].log_mole_frac_phase_comp[Vap,H2]

        Constraints:

            fs.R101.gibbs_minimization[0.0,Vap,CO2]
            fs.R101.control_volume.properties_out[0.0].component_flow_balances[CO2]
            fs.R101.gibbs_minimization[0.0,Vap,CO]
            fs.R101.gibbs_minimization[0.0,Vap,CH3OH]
            fs.R101.control_volume.properties_out[0.0].component_flow_balances[CO]
            fs.R101.control_volume.properties_out[0.0].component_flow_balances[H2O]
            fs.R101.control_volume.properties_out[0.0].component_flow_balances[H2]
            fs.R101.control_volume.properties_out[0.0].log_mole_frac_phase_comp_eqn[Vap,CH3OH]
            fs.R101.control_volume.properties_out[0.0].total_flow_balance
            fs.R101.control_volume.properties_out[0.0].log_mole_frac_phase_comp_eqn[Vap,CO]
            fs.R101.gibbs_minimization[0.0,Vap,H2O]
            fs.R101.control_volume.properties_out[0.0].equilibrium_constraint[Vap,Liq,CH3OH]
            fs.R101.control_volume.properties_out[0.0].component_flow_balances[CH3OH]
            fs.R101.control_volume.properties_out[0.0].sum_mole_frac
            fs.R101.control_volume.properties_out[0.0].log_mole_frac_phase_comp_eqn[Vap,H2O]
            fs.R101.gibbs_minimization[0.0,Liq,CH3OH]
            fs.R101.control_volume.properties_out[0.0]._teq_constraint_Vap_Liq
            fs.R101.control_volume.properties_out[0.0].equilibrium_constraint[Vap,Liq,H2O]
            fs.R101.control_volume.properties_out[0.0].eq_temperature_dew[Vap,Liq,H2O]
            fs.R101.gibbs_minimization[0.0,Liq,H2O]
            fs.R101.control_volume.properties_out[0.0].eq_temperature_dew[Vap,Liq,CO]
            fs.R101.control_volume.properties_out[0.0].eq_temperature_dew[Vap,Liq,CO2]
            fs.R101.control_volume.properties_out[0.0].log_mole_frac_tdew_eqn[Vap,Liq,H2O]
            fs.R101.control_volume.properties_out[0.0].log_mole_frac_comp_eqn[H2O]
            fs.R101.control_volume.properties_out[0.0].eq_mole_frac_tdew[Vap,Liq]
            fs.R101.control_volume.properties_out[0.0].eq_temperature_dew[Vap,Liq,H2]
            fs.R101.control_volume.properties_out[0.0].log_mole_frac_tdew_eqn[Vap,Liq,CH3OH]
            fs.R101.control_volume.properties_out[0.0].eq_temperature_dew[Vap,Liq,CH3OH]
            fs.R101.control_volume.properties_out[0.0].log_mole_frac_comp_eqn[CH3OH]
            fs.R101.gibbs_minimization[0.0,Vap,H2]
            fs.R101.control_volume.properties_out[0.0].log_mole_frac_phase_comp_eqn[Vap,CO2]
            fs.R101.control_volume.properties_out[0.0].log_mole_frac_phase_comp_eqn[Vap,H2]
            fs.R101.control_volume.properties_out[0.0].log_mole_frac_phase_comp_eqn[Liq,H2O]
            fs.R101.control_volume.properties_out[0.0].log_mole_frac_phase_comp_eqn[Liq,CH3OH]

    Independent Block 1:

        Variables:


        Constraints:

            fs.R101.control_volume.properties_out[0.0].sum_mole_frac_out

====================================================================================

[andrewlee94]

The results of the Dulmage-Mendehlson analysis would be a good place to start. Note the last part - Independent Block 1 contains no variable and one constraint (sum of mole fractions), which indicates that the problem is structurally singular; specifically, all the outlet mole fractions are determined by another constraint in the model.

I think I know the issue however - your property package contains constraints for vapor-liquid equilibrium which cannot be combined with the Gibbs energy minimisation (basically, the Gibbs energy minimisation will calculate the phase equilibrium for you). Try removing or deactivating the phase equilibrium constraint (fs.R101.control_volume.properties_out[0.0].equilibrium_constraint) and see what happens.

[panyaoyu]

The results of the Dulmage-Mendehlson analysis would be a good place to start. Note the last part - Independent Block 1 contains no variable and one constraint (sum of mole fractions), which indicates that the problem is structurally singular; specifically, all the outlet mole fractions are determined by another constraint in the model.

I think I know the issue however - your property package contains constraints for vapor-liquid equilibrium which cannot be combined with the Gibbs energy minimisation (basically, the Gibbs energy minimisation will calculate the phase equilibrium for you). Try removing or deactivating the phase equilibrium constraint (fs.R101.control_volume.properties_out[0.0].equilibrium_constraint) and see what happens.

Thank you very much. I am not familiar with GibbsReactor. The reason I used GibbsReactor is that I didn't want to write a reaction package. Now I understand that GibbsReactor might conflict with an inlet flow that contains both gas and liquid phases. I will try to modify the constraints as you suggested, or switch to using YieldReactor.