CalebBell / thermo

Thermodynamics and Phase Equilibrium component of Chemical Engineering Design Library (ChEDL)
MIT License
594 stars 114 forks source link

Need examples for GibbsExcessLiquid + Vapor #104

Open alexchandel opened 2 years ago

alexchandel commented 2 years ago

There are no examples of how to use GibbsExcessLiquid in a flash. This is the most important liquid phase model. An unfinished example is here, but it's unclear from the documentation whether it must be recreated for different temperatures or compositions.

Please provide self-contained examples of flashing with GibbsExcessLiquid, which can be repeatedly evaluated at different temperatures.

alexchandel commented 2 years ago

Here is an example from recent issues:

#!/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+217.0  # this guess is accurate
    pres = 1750000

    res = flasher.flash(zs=zs, P=pres, T=temp)  # use hot_start for speed
    print((res.T, res.H()))

if __name__ == '__main__':
    example()
CalebBell commented 2 years ago

Hi Alex (and anyone else reading this):

The GibbsExcessLiquid class has a few kinks to work out and feature to implement; it is for this reason I didn't create an example.

The only issue in your example above is that the tau parameters are supposed to be either all lists-of-lists, or numpy arrays. The use of tuples is not supported; although they kind of work with the current code, they may not in the future (and some things, like the repr() of the object, do not.

I want to implement Henry's law, and as well quite a few derivatives are missing for certain setting combinations.

Sincerely, Caleb