Help with calculating LLEs

Dear Gustavo,

I hope you are doing great! I am trying to use your really cool package phasepy to calculate some LLEs. I have parameters for both NRTL and UNIQUAC. However, I am not able to reproduce the LLE data points from either of those models using phasepy.

For instance, I have the following info from a binary LLE:

Compound 1: METHANE, TETRACHLORO
- r1 = 3.39
- q1 = 2.91
Compound 2: CYCLOHEXANE, METHYL, PERFLUORO
- r2 = 7.0735
- q2 = 6.44
a12_uniquac = 37.21 @ T=278.15
a21_uniquac = 155.16 @ T=278.15
T: 278.15 K
At this conditions, the LLE compositions should be:
- x1 in phase 1 = 0.295
- x1 in phase 2 = 0.9478

Notice that the expected equilibrium compositions should match exactly the prediction from UNIQUAC because they are smoothed values from UNIQUAC predictions, and correspond to a single-solution Gibbs energy of mixing curve. I got all the data from Vol. 5 of the DECHEMA Chemistry data series. I am able to reproduce the expected equilibrium compositions with my inhouse implementation of UNIQUAC and by solving the isoactivity criterion. However, my implementation depends a lot on the initial guesses of the compositions. This is the reason why I am trying phasepy as I discovered that many of these things have been coded already by you.

I copy my current code here:

from phasepy import component, mixture, virialgamma
from phasepy.equilibrium import lle, lle_init

methane_tetrachloro = component(name='methane_tetrachloro', ri=3.39, qi=2.91)
cyclohexane_methyl_perfluoro = component(name='cyclohexane_methyl_perfluoro', ri=7.0735, qi=6.44)
mix = mixture(methane_tetrachloro, cyclohexane_methyl_perfluoro)
a0 = np.array([[0, 37.21],
              [155.16, 0]])
a1 = np.array([[0, 0],
             [0, 0]])
mix.uniquac(a0, a1)
phase_equilibrium_model = virialgamma(mix, virialmodel='ideal_gas', actmodel='uniquac')
T = 278.15 # K
P = 1.01325 # bar
x0_ref = np.array([0.29, .9])
x0_I, x0_II = lle_init(x0_ref, T, P, phase_equilibrium_model)
LLE = lle(x0_I, x0_II, x0_ref, T, P, phase_equilibrium_model, full_output=True)
x_I_eq, x_II_eq = LLE.X

Would you mind taking a look at my code and see whether I am doing something wrong? I hope you find sometime to help me with this problem...I am super happy to contribute back to phasepy whatever I implement.Thank you so much Gustavo!

Dear Edgar,

I'm sorry for not getting back to you sooner. I've been away the last week. I've been looking at your code, and it seems correct. I do have some questions, though:

Can you check the units of the energy parameters? In phasepy, for the UNIQUAC model, A0 should have units of [K], and A1 should be dimensionless. I also noticed that your global composition ('x0_ref') doesn't add to 1. I suggest sticking to molar fractions that add to 1. It might also be the case that your global composition is close to the outside of the phase boundary. For that reason, the solver might converge to a single phase result. For example, if you use instead the following global composition:

z = np.array([0.6, 0.4]) x0, w0 = lle_init(z, T, P, phase_equilibrium_model) lle(x0, w0, z, T, P, phase_equilibrium_model, full_output=True)

In this case, the results look like this:

       T: 278.15
       P: 1.01325
       error_outer: 1.300357064547822e-09
       error_inner: 1.1451828020774934e-10
       iter: 11
       beta: array([0.46718629, 0.53281371])
       tetha: array([0.])
       X: array([[0.94779661, 0.05220339],
       [0.295042  , 0.704958  ]])
       v: [None, None]
       states: ['L', 'L']

In this case, you have the two liquid phases, and the results are close to the experimental data you mentioned.

I hope that helps.

Kind regards, Gustavo

gustavochm / phasepy

Help with calculating LLEs #20