[SAFT-gamma mie] Errors with LLE calculation

[YouHyin]

Hello, Gustavo!
Your code has been helpful for my research, Thank you. However, I have encountered two major errors during LLE calculations:

In sgtpy\gammamie_mixtures\ares.py:491, while using Numpy's solve function to compute Xass, I encounter a "Singular matrix" error.

import numpy as np
from sgtpy import component, mixture, saftgammamie
from sgtpy.equilibrium import lle, tpd_minimas

def calculate_w_lle(amine_GC):
    # Define components using group contributions
    amine = amine_GC
    water = component(GC={'H2O': 1})
    
    # Create mixture and setup the SAFT-γ-Mie equation of state
    mix = amine + water
    mix.saftgammamie()
    eos = saftgammamie(mix)
    
    # Define molecular weights for mass fraction calculations (g/mol)
    MW_amine = eos.Mw[0]  # Molecular weight of TEA
    MW_H2O = eos.Mw[1]  # Molecular weight of water
    
    # Operating conditions
    T = 350.0  # temperature in K
    P = 1.01325e5  # pressure in Pa
    
    # Set up a range for z[0] values
    z0_values = np.linspace(1e-4, 1-1e-4, 1000)
    
    # Loop through z0_values to find the point where LLE starts to form
    for z0 in z0_values:
        z = np.array([z0, 1-z0])
        x0s, tpd0s = tpd_minimas(2, z, T, P, eos, stateW='L', stateZ='L')
        x0 = x0s[0]
        w0 = x0s[1]
        
        sol = lle(x0, w0, z, T, P, eos, full_output=True)
        x1, x2 = sol.X
        diff = np.abs(x1[0] - x2[0])
        
        if diff > 1e-3:
            # Calculate W_TEA at the first point where LLE is detected
            W_lle = z0 * MW_amine / (z0 * MW_amine + (1 - z0) * MW_H2O)
            return W_lle
    
    # If no LLE is formed, return None 
    return 1
    
amine_1=component(GC={'NH2': 1, 'CH2': 2})
W_solvent_optimal = calculate_w_lle(amine_1)

The error messages are as follows:

File ~\anaconda3\lib\site-packages\sgtpy\gammamie_mixtures\ares.py:491, in d2ares_drho(self, x, rho, temp_aux, Xass0)
    488 d2Dijklab_drho = self.kAB_kl * Fklab
    489 d2Dijklab_drho[self.indexABij] *= d2Iijklab_drho[self.indexAB_id]
--> 491 Xass = Xass_solver(rho, xjvk, DIJ, Dijklab, diagasso, Xass)
    492 CIJ = CIJ_matrix(rho, xjvk, Xass, DIJ, Dijklab, diagasso)
    493 dXass = dXass_drho(rho, xjvk, Xass, DIJ, Dijklab, dDijklab_drho, CIJ)

File ~\anaconda3\lib\site-packages\sgtpy\gammamie_mixtures\association_aux.py:152, in Xass_solver(rho, xjvk, DIJ, Dijklab, diagasso, Xass0)
    150 HIJ[diagasso] -= (xjvk + KIJXass)/Xass
    151 for i in range(15):
--> 152     dXass = np.linalg.solve(HIJ, -dQ)
    153     Xnew = Xass + dXass
    155     is_nan = np.isnan(Xnew)

File ~\anaconda3\lib\site-packages\numpy\linalg\linalg.py:409, in solve(a, b)
    407 signature = 'DD->D' if isComplexType(t) else 'dd->d'
    408 extobj = get_linalg_error_extobj(_raise_linalgerror_singular)
--> 409 r = gufunc(a, b, signature=signature, extobj=extobj)
    411 return wrap(r.astype(result_t, copy=False))

File ~\anaconda3\lib\site-packages\numpy\linalg\linalg.py:112, in _raise_linalgerror_singular(err, flag)
    111 def _raise_linalgerror_singular(err, flag):
--> 112     raise LinAlgError("Singular matrix")

LinAlgError: Singular matrix

In gammamie_mixtures\density_solver.py:91, within the density_topliss function, the value of rho_ub is calculated to be smaller than rho_lb, which causes an error.

amine_2=component(GC={'CH3': 1, 'CH2': 4, 'CH': 1, 'C': 1, 'H2O': 1, 'NH': 1}) W_solvent_optimal = calculate_w_lle(amine_2)

The error messages are as follows:

ValueError                                Traceback (most recent call last)
Cell In[9], line 2
      1 amine_2=component(GC={'CH3': 1, 'CH2': 4, 'CH': 1, 'C': 1, 'H2O': 1, 'NH': 1})
----> 2 W_solvent_optimal = calculate_w_lle(amine_2)

Cell In[3], line 33, in calculate_w_lle(amine_GC)
     30 x0 = x0s[0]
     31 w0 = x0s[1]
---> 33 sol = lle(x0, w0, z, T, P, eos, full_output=True)
     34 x1, x2 = sol.X
     35 diff = np.abs(x1[0] - x2[0])

File ~\anaconda3\lib\site-packages\sgtpy\equilibrium\lle.py:79, in lle(x0, w0, Z, T, P, model, v0, Xass0, K_tol, nacc, full_output)
     77     X0 = np.asarray(xes)
     78     beta0 = np.hstack([beta, 0.])
---> 79     out = multiflash(X0, beta0, equilibrio, Z, T, P, model, v, Xass,
     80                      K_tol, nacc, True)
     81     Xm, beta, tetha, v = out.X, out.beta, out.tetha, out.v
     83 X, W = Xm

File ~\anaconda3\lib\site-packages\sgtpy\equilibrium\multiflash.py:190, in multiflash(X0, betatetha, equilibrium, z, T, P, model, v0, Xass0, K_tol, nacc, full_output)
    187     lnphi[i], v[i], Xass[i] = model.logfugef_aux(X[i], temp_aux, P,
    188                                                  state, v[i], Xass[i])
    189     if np.any(np.isnan(lnphi[i])) or np.isnan(v[i]):
--> 190         lnphi[i], v[i], Xass[i] = model.logfugef_aux(X[i], temp_aux, P, state)
    192 lnK = lnphi[0] - lnphi[1:]
    193 error = np.sum((lnK - lnK_old)**2)

File ~\anaconda3\lib\site-packages\sgtpy\gammamie_mixtures\saftgammamie_mixture.py:1299, in saftgammamie_mix.logfugef_aux(self, x, temp_aux, P, state, v0, Xass0)
   1296 else:
   1297     rho0 = 1./v0
-> 1299 rho, Xass = self.density_aux(x, temp_aux, P, state, rho0, Xass0)
   1300 v = 1./rho
   1301 rhomolecular = Na * rho

File ~\anaconda3\lib\site-packages\sgtpy\gammamie_mixtures\saftgammamie_mixture.py:1022, in saftgammamie_mix.density_aux(self, x, temp_aux, P, state, rho0, Xass0)
    994 """
    995 density_aux(x, temp_aux, P, state)
    996 Method that computes the density of the mixture at given composition,
   (...)
   1019     computed fraction of nonbonded sites
   1020 """
   1021 if rho0 is None:
-> 1022     rho, Xass = density_topliss(state, x, temp_aux, P, Xass0, self)
   1023 else:
   1024     rho, Xass = density_newton(rho0, x, temp_aux, P, Xass0, self)

File ~\anaconda3\lib\site-packages\sgtpy\gammamie_mixtures\density_solver.py:91, in density_topliss(state, x, temp_aux, P, Xass0, saft)
     88 bracket = [rho_lb, rho_ub]
     89 if flag == 1:
     90     # Found inflexion point
---> 91     sol_inf = minimize_scalar(dPsaft_fun, args=(x, temp_aux, saft),
     92                               bounds=bracket, method='Bounded',
     93                               options={'xatol': 1e-1})
     94     rho_inf = sol_inf.x
     95     dP_inf = sol_inf.fun

File ~\anaconda3\lib\site-packages\scipy\optimize\_minimize.py:945, in minimize_scalar(fun, bracket, bounds, args, method, tol, options)
    942     if bounds is None:
    943         raise ValueError('The `bounds` parameter is mandatory for '
    944                          'method `bounded`.')
--> 945     res = _minimize_scalar_bounded(fun, bounds, args, **options)
    946 elif meth == 'golden':
    947     res = _recover_from_bracket_error(_minimize_scalar_golden,
    948                                       fun, bracket, args, **options)

File ~\anaconda3\lib\site-packages\scipy\optimize\_optimize.py:2253, in _minimize_scalar_bounded(func, bounds, args, xatol, maxiter, disp, **unknown_options)
   2250 x1, x2 = bounds
   2252 if not (is_finite_scalar(x1) and is_finite_scalar(x2)):
-> 2253     raise ValueError("Optimization bounds must be finite scalars.")
   2255 if x1 > x2:
   2256     raise ValueError("The lower bound exceeds the upper bound.")

ValueError: Optimization bounds must be finite scalars.

Could you provide some advice on how to resolve these issues?
Thank as always for your help.
Best wish!

[gustavochm]

Hi YouHyin,

I'm very sorry. I received the email about this issue while traveling, and then he forgot about it.
I'll give it a look and see what we can do.

Gustavo

[gustavochm]

Hi YouHyin,

I have been looking into this and there is a problem when the initial guesses are not that great with the lle function. This function has some stability variables that caused some overflow during the LLE calculation. In an updated version of the code I'll add an option to cap the values of those variables, as we only care if they are zero or positive.

For your problem, I'd do the following:

1. First get initial guesses from the tpd minimization
2. Perform the LLE using the flash function. The advantage of this is that in this flash you don't care about the stability vairbales.
3. If the flash function detects LLE, then use the lle function, just to verify the two liquids phases are stable.
4. Only report the LLE if the two liquids phases are different and stable.

Here is the updated function:

def calculate_w_lle(amine_GC):
    # Define components using group contributions
    amine = amine_GC
    water = component(GC={'H2O': 1})
    
    # Create mixture and setup the SAFT-γ-Mie equation of state
    mix = amine + water
    mix.saftgammamie()
    eos = saftgammamie(mix)
    
    # Define molecular weights for mass fraction calculations (g/mol)
    MW_amine = eos.Mw[0]  # Molecular weight of TEA
    MW_H2O = eos.Mw[1]  # Molecular weight of water
    
    # Operating conditions
    T = 350.0  # temperature in K
    P = 1.01325e5  # pressure in Pa
    
    # Set up a range for z[0] values
    z0_values = np.linspace(1e-4, 1-1e-4, 30)
    
    # Loop through z0_values to find the point where LLE starts to form
    for z0 in z0_values:
        z = np.array([z0, 1-z0])
        x0s, tpd0s = tpd_minimas(2, z, T, P, eos, stateW='L', stateZ='L')
        x0 = x0s[0]
        w0 = x0s[1]
        
        # perform flash only is x0 != w0
        if not np.allclose(x0, w0, rtol=1e-3):
            sol = flash(x0, w0, 'LL', z, T, P, eos, full_output=True)
            x1 = sol.X
            x2 = sol.Y
            
            # perform lle only if x1 != x2
            if not np.allclose(x1, x2, rtol=1e-3):
                sol = lle(x1, w2, z, T, P, eos, full_output=True)
                x1, x2 = sol.X
                beta = sol.beta
                tetha = sol.beta
                
                # checking if the result is actually in LLE
                # beta -> phase fraction (ideally bigger than zero for all the phases)
                # tetha -> stability variable (zero for stable phases)
                is_lle = not np.allclose(x1, x2, rtol=1e-3) and np.all(beta > 0) and np.all(tetha == 0.)
        
                if is_lle:
                    # Calculate W_TEA at the first point where LLE is detected
                    W_lle = z0 * MW_amine / (z0 * MW_amine + (1 - z0) * MW_H2O)
                    return W_lle
            
    # If no LLE is formed, return None 
    return 1

Then, this following two codes returns 1 (no LLE). The with np.errstate(all='ignore') is just to avoid all the warnings.

amine_1=component(GC={'NH2': 1, 'CH2': 2})
with np.errstate(all='ignore'):
    W_solvent_optimal = calculate_w_lle(amine_1)
W_solvent_optimal

and

amine_2=component(GC={'CH3': 1, 'CH2': 4, 'CH': 1, 'C': 1, 'H2O': 1, 'NH': 1})
with np.errstate(all='ignore'):
    W_solvent_optimal = calculate_w_lle(amine_2)
W_solvent_optimal

Let me know if that helps,
Gustavo