Question about degrees of freedom

[KevinChengKaiLai]

Hello,

I want to know if there is any method we can use to see what we should specified for degrees of freedom?
for example, we need to set a conversion constraint, how ever we can fix other variable to let DoF = 0.
I was struggling in this specification question for self defined ammonia production process.

Thanks!!

[dallan-keylogic]

There's no general method to determine what needs to be specified to achieve zero degrees of freedom. Generally, if the state variables on inlet ports and some set of performance variables on the unit model are fixed, the unit model should have zero degrees of freedom. (If the inlet port is connected by an Arc to some other port, you should explicitly fix those state variables only during initialization and debugging.)

I think you're using the StoichiometricReactor. On the documentation page, it says that there is one degree of freedom per reaction in the reaction package (blk.rate_reaction_extent) plus the external heat duty when the unit model is configured with has_heat_transfer=True (blk.heat_duty). Therefore, you need to either fix those variables or create constraints to calculate them.

For unit model initialization, the initialization routine might expect certain variables to be fixed. You should examine the initialization routine you're using to determine which, if any, variables are expected to be fixed.

While there may not be one single set of "right" variables to fix, sometimes we end up fixing a wrong set of variables. The Diagnostics Toolbox can help determine if a model is structurally singular via the report_structural_issues method and if a model is numerically singular via the report_numerical_issues method. The SVD Toolbox can help you determine why a model is numerically singular. You can see these tutorials for details.

[KevinChengKaiLai]

Thank you so much, you opinion is great.

I then found out after I set the thermal config of all my components from
"phase_equilibrium_form": {("Vap", "Liq"): fugacity} to "phase_equilibrium_form": {("Vap"): fugacity},

Also simplify the phases from
"phases": {"Liq": {"type": LiquidPhase, "equation_of_state": Ideal}, "Vap": {"type": VaporPhase, "equation_of_state": Ideal}, }
to
"phases": { "Vap": {"type": VaporPhase, "equation_of_state": Ideal}, },

The degree of freedom problem is solved, so I suppose it is the liquid-vapor equilibrium increasing DOF.

I originally fix the feed as below (I am trying to design an simple ammonia production)
m.fs.Feed.outlet.flow_mol_phase_comp[0, "Vap", "nitrogen"].fix( 250 * pyunits.mol / pyunits.s ) m.fs.Feed.outlet.flow_mol_phase_comp[0, "Vap", "hydrogen"].fix( 750 * pyunits.mol / pyunits.s ) m.fs.Feed.outlet.flow_mol_phase_comp[0, "Vap", "ammonia"].fix( 1e-5 * pyunits.mol / pyunits.s ) m.fs.Feed.outlet.flow_mol_phase_comp[0, "Liq", "nitrogen"].fix( 1e-5 * pyunits.mol / pyunits.s ) m.fs.Feed.outlet.flow_mol_phase_comp[0, "Liq", "hydrogen"].fix( 1e-5 * pyunits.mol / pyunits.s ) m.fs.Feed.outlet.flow_mol_phase_comp[0, "Liq", "ammonia"].fix( 1e-5 * pyunits.mol / pyunits.s )

and the DOF was still high, can you provide any suggestion where I miss to fix if I want to design both Liquid and Vapor reaction system ?

[dallan-keylogic]

You aren't including enough information for me to fully troubleshoot your model. I would need to see the full configuration you're passing to the modular property framework as well as the reaction package you're using for ammonia synthesis. Also, I don't know what units are present in your system.

First, it looks like you're using FpcTX state variables. That set of state variables has issues with degenerate equations when there is only a single phase present. I'd recommend switching to FTPx for most of your flowsheet, perhaps only using FpcTX for the flash vessel itself (assuming you're trying to set up a conventional Haber-Bosch style process) and converting between that and FTPx via TranslatorBlocks.

Second, hydrogen and nitrogen shouldn't be liquid phase components (except in extraordinary cases where you're concerned about tiny concentrations of dissolved hydrogen). Look at this implementation of Peng Robinson for natural gas components to see how the valid_phase_types argument works.

Third, "phase_equilibrium_form": {("Vap"): fugacity} is incorrect. You specify a form for phase equilibrium pairs, like ("Vap", "Liq"). It might not raise an error because the property package doesn't bother looking for phase equilibrium forms when there is only a single phase present.

Finally, have you gone through your flowsheet, fixed the variables on your inlet ports, and checked whether the resulting unit models have zero degrees of freedom? That process might help you find which unit model has the DOF issue.

[KevinChengKaiLai]

Hi, sorry for the inconvenience,

For now, I didn't know how to use FTPx, so I'm still using FpcTP in my configuration (as the idaes example does).
I did make the degrees of freedom become 0 but I think that is in the wrong situation since my configuration file is not correct so I will try to redo it.
below is the thermal configuration:

thermal_config = {
    # Specifying components
    "components": {
        "nitrogen": {
            "type": Component,
            "elemental_composition": {"N": 2},
            "dens_mol_liq_comp": Perrys,
            "enth_mol_liq_comp": Perrys,
            "enth_mol_ig_comp": RPP4,
            "pressure_sat_comp": RPP4,
            "phase_equilibrium_form": {("Vap", "Liq"): fugacity},
            "parameter_data": {
                "mw": (28.0134-3, pyunits.kg / pyunits.mol),  # [1] pg. 676
                "pressure_crit": (3.39e6, pyunits.Pa),  # [1] pg. 676
                "temperature_crit": (126.2, pyunits.K),  # [1] pg. 676
                "dens_mol_liq_comp_coeff": {  # [2] pg. 2-97 
                    "eqn_type": 1,
                    "1": (3.2091, pyunits.kmol * pyunits.m**-3),
                    "2": (0.2861, None),
                    "3": (126.2, pyunits.K),
                    "4": (0.2966, None),
                },
                # DONE
                "cp_mol_ig_comp_coeff": {  # [1] pg. 664
                    "A": (3.115e1, pyunits.J / pyunits.mol / pyunits.K),
                    "B": (-1.357e-2, pyunits.J / pyunits.mol / pyunits.K**2),
                    "C": (2.680e-5, pyunits.J / pyunits.mol / pyunits.K**3),
                    "D": (-1.168e-8, pyunits.J / pyunits.mol / pyunits.K**4),
                },
                # DONE
                "cp_mol_liq_comp_coeff": {  # [2] pg. 2-173
                    "1": (2.8197e5, pyunits.J / pyunits.kmol / pyunits.K),
                    "2": (-1.2281e4, pyunits.J / pyunits.kmol / pyunits.K**2),
                    "3": (2.4800e2, pyunits.J / pyunits.kmol / pyunits.K**3),
                    "4": (-2.2182, pyunits.J / pyunits.kmol / pyunits.K**4),
                    "5": (7.4902e-3, pyunits.J / pyunits.kmol / pyunits.K**5),
                },
                # DONE
                "enth_mol_form_liq_comp_ref": (
                    8.669e3,
                    pyunits.J / pyunits.mol,
                ),  # [3] updated 5/10/24
                # DONE
                "enth_mol_form_vap_comp_ref": (
                    0, # element delta H = 0
                    pyunits.J / pyunits.mol,
                ),  # [3] updated 5/10/24
                # DONE
                "pressure_sat_comp_coeff": {
                    "A": (-6.09676, None),  # [1] pg. 666
                    "B": (1.13670, None),
                    "C": (-1.04072, None),
                    "D": (-1.93306, None),
                },
            },
        },
        "hydrogen": {
            "type": Component,
            "elemental_composition": {"H": 2},
            "dens_mol_liq_comp": Perrys,
            "enth_mol_liq_comp": Perrys,
            "enth_mol_ig_comp": RPP4,
            "pressure_sat_comp": RPP4,
            "phase_equilibrium_form": {("Vap", "Liq"): fugacity},
            "parameter_data": {
                "mw": (2.016e-3, pyunits.kg / pyunits.mol),  # [1] pg. 667
                "pressure_crit": (13e6, pyunits.Pa),  # [1] pg. 667
                "temperature_crit": (32.97, pyunits.K),  # [1] pg. 667
                "dens_mol_liq_comp_coeff": {  # [2] pg. 2-97 
                    "eqn_type": 1,
                    "1": (5.414, pyunits.kmol * pyunits.m**-3),
                    "2": (0.34893, None),
                    "3": (33.19, pyunits.K),
                    "4": (0.2706, None),    
                },
                # DONE
                "cp_mol_ig_comp_coeff": {  # [1] pg. 668 CPVAP A B C D
                    "A": (2.714e1, pyunits.J / pyunits.mol / pyunits.K),
                    "B": (9.274e-3, pyunits.J / pyunits.mol / pyunits.K**2),
                    "C": (-1.381e-5, pyunits.J / pyunits.mol / pyunits.K**3),
                    "D": (7.645e-9, pyunits.J / pyunits.mol / pyunits.K**4),
                },
                # DONE
                "cp_mol_liq_comp_coeff": {  # [2] pg. 2-174
                    "1": (6.6653e1, pyunits.J / pyunits.kmol / pyunits.K),
                    "2": (6.7659e3, pyunits.J / pyunits.kmol / pyunits.K**2),
                    "3": (-1.2363e2, pyunits.J / pyunits.kmol / pyunits.K**3),
                    "4": (4.7827e2, pyunits.J / pyunits.kmol / pyunits.K**4),
                    "5": (0, pyunits.J / pyunits.kmol / pyunits.K**5),
                },
                # DONE -- Question how to find liquid molar enthalpy
                "enth_mol_form_liq_comp_ref": (
                    7.9252e3,
                    pyunits.J / pyunits.mol,
                ),  # [3] updated 5/10/24
                # DONe
                "enth_mol_form_vap_comp_ref": (
                    0,
                    pyunits.J / pyunits.mol,
                ),  # [3] updated 5/10/24
                "pressure_sat_comp_coeff": {
                    "A": (-5.57929, None),  # [1] pg. 664 VP A B C D
                    "B": (2.60012, None),
                    "C": (-0.85506, None),
                    "D": (1.70503, None),
                },
            },
        },
        "ammonia": {
            "type": Component,
            "elemental_composition": {"H": 3, "N": 1},
            "dens_mol_liq_comp": Perrys,  # fitted to this equation form
            "enth_mol_liq_comp": Perrys,  # fitted to this equation form
            "enth_mol_ig_comp": NIST,
            "pressure_sat_comp": RPP4,  # fitted to this equation form
            "phase_equilibrium_form": {("Vap", "Liq"): fugacity},
            "parameter_data": {
                "mw": (17.031e-3, pyunits.kg / pyunits.mol),  # [4] 5/20/24
                "pressure_crit": (11.3e6, pyunits.Pa),  # [4] 5/20/2024
                "temperature_crit": (405.56, pyunits.K),  # [4] 5/202/24
                "dens_mol_liq_comp_coeff": {  # [2] pg. 2-97 DONE
                    "eqn_type": 1,
                    "1": (3.5383, pyunits.kmol * pyunits.m**-3),
                    "2": (0.25443, None),
                    "3": (405.65, pyunits.K),
                    "4": (0.2888, None),    
                },
                "cp_mol_ig_comp_coeff": {  # DONE [3] valid on 298-1200 K, updated 5/10/2024
                    "A": (19.99563, pyunits.J / pyunits.mol / pyunits.K),
                    "B": (49.77119, pyunits.J / pyunits.mol / pyunits.K / pyunits.kK),
                    "C": (-15.37599,
                        pyunits.J / pyunits.mol / pyunits.K / pyunits.kK**2,
                    ),
                    "D": (
                        1.921168,
                        pyunits.J / pyunits.mol / pyunits.K / pyunits.kK**3,
                    ),
                    "E": (
                        0.189174,
                        pyunits.J / pyunits.mol / pyunits.K / pyunits.kK**-2,
                    ),
                    "F": (-53.30667, pyunits.kJ / pyunits.mol),
                    "G": (203.8591, pyunits.J / pyunits.mol / pyunits.K),
                    "H": (-45.89806, pyunits.kJ / pyunits.mol),
                },
                # DONE
                "cp_mol_liq_comp_coeff": {  # [6] pg. 1189 regressed from x1=1 and Cp/R for all T values
                    "1": (6.1289e1, pyunits.J / pyunits.kmol / pyunits.K),
                    "2": (8.0925e4, pyunits.J / pyunits.kmol / pyunits.K**2),
                    "3": (7.994e2, pyunits.J / pyunits.kmol / pyunits.K**3),
                    "4": (-2.6510e3, pyunits.J / pyunits.kmol / pyunits.K**4),
                    "5": (0, pyunits.J / pyunits.kmol / pyunits.K**5),
                },
                # DONE -- perplexity pro
                "enth_mol_form_liq_comp_ref": (
                    -80.882e3,
                    pyunits.J / pyunits.mol,
                ),  # [6] Table 4 pg. 1183
                "enth_mol_form_vap_comp_ref": (
                    -45.9e3,
                    pyunits.J / pyunits.mol,
                ),  # [3] updated 5/10/24
                "pressure_sat_comp_coeff": {
                    "A": (-18.122, None),  # TODO [5] data for regression on pg. 6-122
                    "B": (10.596, None),
                    "C": (-18.908, None),
                    "D": (-32.728, None),
                },
            },
        },
        
    },
    # Specifying phases
    "phases": {
        # "Liq": {"type": LiquidPhase, "equation_of_state": Ideal},
        "Vap": {"type": VaporPhase, "equation_of_state": Ideal},
    },
    # Set base units of measurement
    "base_units": {
        "time": pyunits.s,
        "length": pyunits.m,
        "mass": pyunits.kg,
        "amount": pyunits.mol,
        "temperature": pyunits.K,
    },
    # Specifying state definition
    "state_definition": FpcTP,
    "state_bounds": {
        "flow_mol_phase_comp": (0, 100, 1000, pyunits.mol / pyunits.s),
        "temperature": (273.15, 298.15, 1000, pyunits.K),
        "pressure": (5e4, 1e5, 1e7, pyunits.Pa),
    },
    "pressure_ref": (1e5, pyunits.Pa),
    "temperature_ref": (298.15, pyunits.K),
    # Defining phase equilibria
    # "phases_in_equilibrium": [("Vap", "Liq")],
    # "phase_equilibrium_state": {("Vap", "Liq"): SmoothVLE},
    # "bubble_dew_method": IdealBubbleDew,
    # "default_scaling_factors": {"mole_frac_comp[benzene]": 1},
}

below is the reaction configuration

reaction_config = {
    "base_units": {
        "time": pyunits.s,
        "length": pyunits.m,
        "mass": pyunits.kg,
        "amount": pyunits.mol,
        "temperature": pyunits.K,
    },
    "rate_reactions": {
        "NH3_synthesis": {
            "stoichiometry": {
                ("Vap", "nitrogen"): -1,
                ("Vap", "hydrogen"): -3,
                ("Vap", "ammonia"): 2
            },
            "heat_of_reaction": constant_dh_rxn,
            "rate_constant": arrhenius,
            "rate_form": power_law_rate,
            "concentration_form": ConcentrationForm.molarity,  # m3 or L?
            "parameter_data":{

                "reaction_order": {("Vap", "hydrogen"): 1,
                                ("Vap", "nitrogen"): 1}, 
                "dh_rxn_ref": (-92.28e3, pyunits.J / pyunits.mol), # DONE
                "arrhenius_const": (1.5638e-9, 
                                    pyunits.m**3 / pyunits.mol * pyunits.s**-1), # TODO
                "energy_activation": (40,765, pyunits.J / pyunits.mol), # DONE
            }
        }    
    }
}