Equilibrium calculations

Different equilibrium calculations could be performed with the yaeos Excess Gibbs models.

Let’s instantiate an example model.

Example model

As always, all the described methods are available for all Excess Gibbs models. But we are going to use an ethanol-hexane-water mixture modeled with UNIFAC as example:

import yaeos


functional_groups = [
    {1: 1, 2: 1, 14: 1},    # ethanol groups
    {1: 2, 2: 4},           # n-hexane groups
    {16: 1},                # water group
]

model = yaeos.UNIFACVLE(functional_groups)

Flash \(\gamma \text{-} \gamma\)

\[y_{i} \gamma_{i}^{y} = x_{i} \gamma_{i}^{x}\]

Liquid-liquid flash calculation with the \(\gamma \text{-} \gamma\) approach could be performed as follows:

import numpy as np


T = 30 + 273.15             # Temperature [K]
n = np.array([5, 10, 10])   # Moles of each substance [mol]

# IMPORTANT: flash calculations are performed specifying mole fractions
z = n / n.sum()             # Mole fraction of eac substance [-]


flash = model.flash_t(z, T)

flash

{'x': array([0.03080572, 0.96735068, 0.0018436 ]),
 'y': array([0.31604605, 0.01086864, 0.67308532]),
 'T': 303.15,
 'beta': 0.5931639638170902}

The result of the flash calculation is returned as a Python dictionary with the keys:

• x: mole fractions of phase \(x\)

• y: mole fractions of phase \(y\)

• T: temperature in Kelvin that was specified

• beta: mole fraction of phase \(y\)

We can calculate the activities on each phase to check that are equal:

x = flash["x"]
y = flash["y"]
beta = flash["beta"]

ln_gammas_x = model.ln_gamma(x, T)
ln_gammas_y = model.ln_gamma(y, T)

gammas_x = np.exp(ln_gammas_x)
gammas_y = np.exp(ln_gammas_y)

activities_x = gammas_x * x
activities_y = gammas_y * y

print(activities_x)
print(activities_y)

[0.46710663 0.97593751 0.87010966]
[0.46710662 0.9759375  0.87010961]

As shown the activities in both phases are equals, the flash calculation was satisfactory.

We can use the beta value to obtain the moles of each phase:

nty = beta * n.sum()         # Total moles of phase y
ntx = (1 - beta) * n.sum()   # total moles of phase x

ny = y * nty  # moles vector of phase y
nx = x * ntx  # moles vector of phase x

print("Mixutre: ethanol-hexane-water")
print("Global mole vector: ", n)
print("Phase y mole vector: ", ny)
print("Phase x mole vector: ", nx)

Mixutre: ethanol-hexane-water
Global mole vector:  [ 5 10 10]
Phase y mole vector:  [4.68667812 0.16117207 9.9812489 ]
Phase x mole vector:  [0.31332188 9.83882793 0.0187511 ]

There is an additional parameter to the flash_t, k0, which allows to provide an initialization to the flash calculation with:

\[k_{0} = \frac{y}{x}\]

The complete documentation of the flash_t can be accessed with:

help(model.flash_t)

Help on method flash_t in module yaeos.core:

flash_t(z, temperature: float, k0=None) -> dict method of yaeos.models.excess_gibbs.unifac.UNIFACVLE instance
    Two-phase split with specification of temperature and pressure.

    Parameters
    ----------
    z : array_like
        Global mole fractions
    temperature : float
        Temperature [K]
    k0 : array_like, optional
        Initial guess for the split, by default None (will use stability
        analysis)

    Returns
    -------
    dict
        Flash result dictionary with keys:
            - x: "phase 1" mole fractions
            - y: "phase 2" mole fractions
            - T: temperature [K]
            - beta: "phase 2" fraction

Stability analysis

It’s possible to perform stability analysis by tm minimization proposed by Michelsen.

Let’s check with an example:

T = 30 + 273.15             # Temperature [K]
n = np.array([5, 10, 10])   # Moles vector [mol]

z = n / n.sum()  # global mole fractions


tm_min, tm_all_min = model.stability_analysis(z, T)

print(tm_min)

{'w': array([2.07328795e-02, 1.47170347e-04, 9.79119950e-01]), 'tm': -0.5172683567506908}

print(tm_all_min)

{'tm': array([-0.51726836, -0.38462711, -0.51726836]), 'w': array([[2.07328795e-02, 1.47170347e-04, 9.79119950e-01],
       [6.54275506e-03, 9.92256323e-01, 1.20092189e-03],
       [2.07328782e-02, 1.47170324e-04, 9.79119951e-01]])}

Three minimization were performed, starting from near pure compositions of each compound.

The first output tm_min, provide the value and composition of the lowest of the all found minimums.

If the tm value of tm_min is negative, the phase of composition z is considered unstable. If you check the used composition, is the same used in the flash \(\gamma \text{-} \gamma\) tutorial. There, was shown that a flash calculation converges under the same conditions. The stability analysis in this case:

tm_min["tm"]

-0.5172683567506908

has a negative value, confirming that the system was truly unstable.

In the second output tm_all_min are listed all the found minimums. You can check that the tm_min is also listed there and repeated. Different initialization could converge to the same minimum.

You can check in the minimums list that there is another minimum with value -0.38462711. We can use the composition of the two found minimums as initialization of a flash calculation:

y = tm_all_min["w"][0]
x = tm_all_min["w"][1]

k0 = y / x

flash = model.flash_t(z, T, k0)

flash

{'x': array([0.03080571, 0.96735068, 0.0018436 ]),
 'y': array([0.31604605, 0.01086864, 0.67308532]),
 'T': 303.15,
 'beta': 0.5931639641540031}