Flashing with H2 throws nonsense `ValueError: math domain error`

[alexchandel]

Flashing at intermediate temperatures with hydrogen (e.g. 600 C) results in the following nonsense error:

./thermo_error.py
Traceback (most recent call last):
  File "/Users/alex/Desktop/./thermo_error.py", line 90, in <module>
    example()
  File "/Users/alex/Desktop/./thermo_error.py", line 86, in example
    print(flasher.flash(zs=zs, P=pres, T=temp))
  File "/usr/local/lib/python3.9/site-packages/thermo/flash/flash_base.py", line 275, in flash
    g, ls, ss, betas, flash_convergence = self.flash_TPV(T=T, P=P, V=V, zs=zs, solution=solution, hot_start=hot_start)
  File "/usr/local/lib/python3.9/site-packages/thermo/flash/flash_vl.py", line 746, in flash_TPV
    return self.flash_TP_stability_test(T, P, zs, self.liquid, self.gas, solution=solution)
  File "/usr/local/lib/python3.9/site-packages/thermo/flash/flash_vl.py", line 640, in flash_TP_stability_test
    stable, (trial_zs, appearing_zs, V_over_F, stab_guess_name, stab_guess_number, stab_sum_zs_test, stab_lnK_2_tot) = self.stability_test_Michelsen(
  File "/usr/local/lib/python3.9/site-packages/thermo/flash/flash_vl.py", line 549, in stability_test_Michelsen
    lnK = log(Ks[k])
ValueError: math domain error

Self-contained example code below replicates this error. Removing the hydrogen eliminates the error. Other gas models (e.g. IdealGas) have no effect; the issue appears to be that hydrogen in the (nonexistent) liquid phase calculation.

Even though a human could tell no liquid phase exists at this temperature, this issue must still be fixed due to dependent issues like #106, where thermo's solvers test such temperatures. And generally, even if no liquid phase exists, thermo should not crash.

#!/usr/bin/env python

from thermo import ChemicalConstantsPackage, PRMIX, CEOSGas, FlashVL, \
    EquilibriumState, GibbsExcessLiquid
from thermo.interaction_parameters import IPDB
from thermo.nrtl import NRTL

def get_nrtl(zs: list[float], T: float = 298.15) -> NRTL:
    """Build NRTL."""
    tau_as = (
        (0,         -5.1549,    0,          0,      0),
        (5.8547,    0,          -0.40926,   0,      0),
        (0,         -0.39036,   0,          0,      0),
        (0,         0,          0,          0,      0),
        (0,         0,          0,          0,      0),
    )
    tau_bs = (
        (0,         2270.62,    284.966,    0,      0),
        (229.497,   0,          1479.46,    0,      0),
        (-216.256,  447.003,    0,          0,      0),
        (0,         0,          0,          0,      0),
        (0,         0,          0,          0,      0),
    )
    alpha_cs = (
        (0,     0.2,    0.3,    0,      0),
        (0.2,   0,      0.46,   0,      0),
        (0.3,   0.46,   0,      0,      0),
        (0,     0,      0,      0,      0),
        (0,     0,      0,      0,      0),
    )
    return NRTL(T=T, xs=zs, tau_as=tau_as, tau_bs=tau_bs, alpha_cs=alpha_cs)

def get_nrtl_flasher(
    names: list[str],
    zs: list[float] = None,
    T: float = 298.15,
    P: float = 101325,
) -> tuple[NRTL, FlashVL]:
    """Build NRTL flasher."""
    if zs is None:
        zs = [1.0/len(names)] * len(names)
    consts, props = ChemicalConstantsPackage.from_IDs(names)
    kijs = IPDB.get_ip_asymmetric_matrix('ChemSep PR', consts.CASs, 'kij')
    eos_kwargs = {'Pcs': consts.Pcs, 'Tcs': consts.Tcs, 'omegas': consts.omegas, 'kijs': kijs}
    gas = CEOSGas(PRMIX, eos_kwargs=eos_kwargs, HeatCapacityGases=props.HeatCapacityGases)

    nrtl = get_nrtl(zs=zs, T=T)
    liquid = GibbsExcessLiquid(
        VaporPressures=props.VaporPressures,
        HeatCapacityGases=props.HeatCapacityGases,
        VolumeLiquids=props.VolumeLiquids,
        GibbsExcessModel=nrtl,
        equilibrium_basis='Psat', caloric_basis='Psat',
        T=T, P=P, zs=zs)

    return FlashVL(consts, props, liquid=liquid, gas=gas)

def add_heat(flasher: FlashVL,
             zs: list[float],
             hrxn: float,
             T0: float,
             P: float = 101325,
             extra_cp: float = 0,) -> EquilibriumState:
    """Flash with Hrxn."""
    eqs: EquilibriumState = flasher.flash(T=T0, P=P, zs=zs)
    current_val = eqs.H()
    cp0 = eqs.dH_dT_P()
    hrxn_vl = hrxn * (cp0 / (cp0 + extra_cp))

    return flasher.flash(zs=zs, P=P, H=current_val + hrxn_vl)

def example():
    """Run."""
    stoich = {'nitrobenzene': -1, 'water': +2, 'aniline': +1, 'hydrogen': -3, 'methane': 0}
    flasher = get_nrtl_flasher(list(stoich.keys()))
    zs = [0, 0.81, 0.16, 0.02, 0.01]
    temp = 273.15+600
    pres = 1750000

    print(flasher.flash(zs=zs, P=pres, T=temp))

if __name__ == '__main__':
    example()

[CalebBell]

Hi Alex If you install thermo==0.2.18 and chemicals==1.0.19 you may have a little more luck, assuming you do not use the hot_start parameter.

[alexchandel]

@CalebBell Thank you! Apparently my env was finding an old thermo==0.2.14 install.