Phase envelope fails for mixtures

[aegis4048]

The below is the reproducible codes. It represents a simple hydrocarbon mixture of 20% methane, 20% propane, and 60% hexane.

` from thermo import ChemicalConstantsPackage, PRMIX, CEOSLiquid, CEOSGas, PropertyCorrelationsPackage from thermo import FlashVLN import matplotlib.pyplot as plt from thermo.interaction_parameters import IPDB

def psi_to_Pa(psi): return psi * 6894.745

def F_to_K(F): return (F - 32) * 5/9 + 273.15

def Pa_to_psi(pa): return pa / 6894.745

def K_to_F(K): return (K - 273.15) * 1.8 + 32

comp = dict([ ('methane', 0.2), ('propane', 0.2), ('hexane', .6), ('water', 0.), ])

total_comp = sum(comp.values()) if total_comp > 1: comp = {k: v / total_comp for k, v in comp.items()}

constants, properties = ChemicalConstantsPackage.from_IDs(comp.keys()) kijs = IPDB.get_ip_asymmetric_matrix('ChemSep PR', constants.CASs, 'kij')

eos_kwargs = dict(Tcs=constants.Tcs, Pcs=constants.Pcs, omegas=constants.omegas, kijs=kijs) gas = CEOSGas(PRMIX, eos_kwargs, HeatCapacityGases=properties.HeatCapacityGases) liq = CEOSLiquid(PRMIX, eos_kwargs, HeatCapacityGases=properties.HeatCapacityGases) flashN = FlashVLN(constants, properties, liquids=[liq, liq], gas=gas) # if dealing with water, fill [liq, liq] twice.

P = psi_to_Pa(100) T = F_to_K(100)

res = flashN.flash(T=T, P=P, zs=list(comp.values()))

zs = list(comp.values())

fig = flashN.plot_TP(zs=zs, show=False) ax = fig.get_axes()[0] ax.set_yscale('linear') `

You can notice that the dew point curve is near-zero. Not only that, the bubble point curve also seems to be wrong. The phase envelope works correctly for pure components, but fails when significant mixture fractions are introduced. This is very weird to me because Thermo correctly calculates vapor fraction and liquid fraction at a given T and P with the flash module... as far as I know if it can correctly predict vapor vs. liquid fraction, it should also be able to easily predict bubble and dew points.

The below is the screenshot of the result I got from the commercial software I'm using for the 0.2 0.2 0.6 mixture:

I found similar issues when plotting isobaric viscosity against temperature (accurate for pure compound, bad for mixtures), but I think they share the same underlying issue.

Is there something that I'm missing here or is this an actual bug? I've been struggling the past 10~15 hrs over this... I realized I can't resolve this on my own so I'm opening an issue here, please help.

[CalebBell]

Hi Eric, In general commercial packages often have quite polished phase envelope plotting modules. It is a difficult task, however. The functionality in Thermo is simple and starts off from the lowest component's melting point, which seems to sometimes solve to a false solution in this case.

I don't have any examples using the function you found, plot_TP, in part because cubic EOSs are difficult to plot in a way activity coefficient modules, or ideal components, or pure components are not.

Thermo can probably plot roughly the same plot as you showed if the same parameters were used, but it's far from perfect. Please remember both the conditions of the license, and my many warnings on using the package. It is a hobby project. Providing an initial guess and using the not-recommended option hot which continues each flash from the previous flash, the plot gets quite close to connecting the curves:

fig = flashN.plot_TP(zs=zs, show=True, Tmin=260, Tmax=475, pts=50, hot=True)

Sincerely, Caleb

[aegis4048]

Thanks for the reply.

I tried setting hot=True and it does seem to improve but as you mentioned it does have some room to improve. I'm trying to implement an improved phase diagram model using some of thermo's available attributes like interaction parameters, will update here when I get it done, planning to make a commit when I'm done if possible.