Calculating and plotting LLEs

[StvnF02]

Hello,
I'm using FeOs to calculate liquid-liquid-equilibria (in python). To get started, I tried reconstructing the VLE/LLE found in 10.1021/ie010954d Fig.8 (methanol - cyclohexane). At first i used feos.eos.PhaseDiagram.binary_vlle, however this returns a PhaseDiagramHetero which apparently does not have a to_dict - function. Is there a different way of plotting a PhaseDiagramHetero other than converting it into a pd.DataFrame and plotting the dataframe? If so, where is it documented?

Afterwards, I used a combination of a feos.eos.PhaseDiagram.binary_vle and a feos.eos.PhaseDiagram.binary_lle, as they both return a PhaseDiagram object which can be converted to a dataframe. When comparing my results to the ones provided in 10.1021/ie010954d Fig.8, the LLE looks off, especially at lower methanol mole fractions, as can be seen in the picture attached (the coloured lines being my results, black lines coming from the paper).

I'm now wondering if the issue originates from a mistake on my end or if perhaps from the LLE calcultation itself? If I did something wrong, please don't hesitate to tell me 😄

Lastly, I'm planning on calculating and plotting ternary LLE in the future. Is that possible with FeOs? I did not find anything about ternary systems (especially plotting them) in the FeOs documentation!

Kind regards

[g-bauer]

Hey,

thanks for the issue. You can access all calculated State objects inPhaseDiagram objects by either using the states or the vapor and liquid fields. Here is an example notebook that explains PhaseDiagram in detail (for a pure substance but the same applies for mixtures).

Here is an example that you can use to generate the plot above (you'll have to adjust the path to the parameter files):

import feos
from feos.pcsaft import PcSaftParameters
from feos.eos import *
import numpy as np
import si_units as si
import matplotlib.pyplot as plt

parameters = PcSaftParameters.from_multiple_json(
    [(["methanol"], "../../parameters_old/pcsaft/gross2002.json"),
    (["cyclohexane"], "../../parameters_old/pcsaft/gross2001.json")],
    "../../parameters_old/pcsaft/gross2002_binary.json"
)
eos = EquationOfState.pcsaft(parameters)
vle = PhaseDiagram.binary_vle(eos, 1.013 * si.BAR, npoints=100)
lle = PhaseDiagram.lle(
    eos, 
    1.013 * si.BAR, 
    feed=np.array([0.5, 0.5])*si.MOL, 
    min_tp=5 * si.CELSIUS, 
    max_tp=50 * si.CELSIUS, npoints=100
)

plt.plot(vle.liquid.molefracs[:, 0], vle.liquid.temperature / si.CELSIUS, color="blue")
plt.plot(vle.vapor.molefracs[:, 0], vle.liquid.temperature / si.CELSIUS, color="blue")

plt.plot(lle.liquid.molefracs[:, 0], lle.liquid.temperature / si.CELSIUS, color="red")
plt.plot(lle.vapor.molefracs[:, 0], lle.liquid.temperature / si.CELSIUS, color="red")

plt.ylabel("T / °C")
plt.xlabel(r"$x_\text{Methanol}$")
plt.ylim(0, 90)
plt.xlim(0, 1);

Here I am using PhaseDiagram.binary_vle for the VLE and PhaseDiagram.lle for the LLE. Note that the states are called vapor and liquid even though they are both liquid phases for the LLE (we might have to rename those).

Once you have a PhaseDiagram object, you can calculate any properties of any state. For example, you can loop over all liquid state points of the VLE:

for state in vle.liquid:
    print(state.molefracs)
    print(state.pressure())
    print(state.molar_enthalpy(Contributions.Residual))
    print(state.molar_entropy(Contributions.Residual))

or, if you need both equilibrium states at once, you can use states to access the PhaseEquilibrium objects (that contain the states of all phases in equilibrium):

for vle in vle.states:
    print(vle.vapor.molar_enthalpy(Contributions.Residual))
    print(vle.liquid.molar_enthalpy(Contributions.Residual))

For now, we don't provide any plotting functionalities but since you can access states, like shown above, it should be straight forward to create plots for any property FeOs can calculate.