[SAFT-Gamma-Mie] Errors with flash calculation

[YouHyin]

Hello Gustavo!

Your updates to sgtpy have been a great help in my research. Thank you very much. The flash calculations are converging much better than before, but I have a problem.

import numpy as np
import matplotlib.pyplot as plt
from sgtpy import component, mixture, saftgammamie
from sgtpy.equilibrium import flash


# Set the necessary component and mixture settings for the actual calculation
T = 313.15  # K
P = 101000  # Pa

# mass fraction of amine
mass_amine = 0.304

# Initial lean solvent loading
lean_loading=1.5

# Flow rates of lean solvent and inlet gas 
flue_gas_flow = 220.29 # kmol/hr
lean_solvent_flow = 39.20 # kmol/hr


N2 = component(GC={'N2': 1})
CO2 = component(GC={'CO2': 1})
H2O = component(GC={'H2O': 1})
AEEA = component(GC={'NH2': 1, 'CH2': 3, 'NH':1,'CH2OH': 1})

mix = N2 + CO2 + AEEA + H2O
mix.saftgammamie()
eos = saftgammamie(mix)
# For initial lean solvent, set N2 mole fraction to nearly zero
moles_N2 = 1E-6
alpha = lean_loading

mass_amine = mass_amine
moles_amine = mass_amine / eos.Mw[2]

moles_co2 = alpha * moles_amine
mass_co2 = moles_co2 * eos.Mw[1]

mass_base = 1.
mass_water = mass_base - mass_amine - mass_co2

moles_water = mass_water / eos.Mw[3]
moles_array = np.array([moles_N2, moles_co2, moles_amine, moles_water])

# Initial lean amine solvent mole fraction (N2, CO2, amine, water 순서)
x = moles_array / np.sum(moles_array, axis=0)
x0 = x
print(f"x0: {x0}")



#flue gas composition N2, CO2, amine, water
y0 = np.array([0.771,0.119,1E-5,0.110])

# Calculate the mole flows for each component in flue gas and lean solvent
flue_gas_mole_flows = flue_gas_flow*y0
lean_solvent_mole_flows = lean_solvent_flow*x0

# Calculate the total mole flows for each component
total_mole_flows = flue_gas_mole_flows + lean_solvent_mole_flows


# Calculate the global compositions
z = total_mole_flows / np.sum(total_mole_flows)
print(f"z(global compositon): {z}")


sol=flash(x0, y0, 'LV', z, T, P, eos, minimization_approach='helmholtz',not_in_x_list=[0],full_output=True) # Set the amount of N2 in the liquid phase to zero

print(sol.X, sol.Y)

When I run this code, the liquid and vapor mole fractions come out as nan. I think that during the flash calculation process, there might be an formula in the form of 1/x, and x is becoming 0, leading to nan values. I attempted to debug and change this expression to the form 1/(x+10^-8) to prevent it from diverging to infinity, but I couldn't find the exact expression.

Can you give me any advice on how to change the formula in the flash calculation process so that the nan value does not come out?
Thank as always for your help.
Best wish!

[gustavochm]

Dear YouHyin,

Thank you very much for your question! The fact you are testing sgtpy with very challenging mixtures really helps me to keep improving the code. So there are two things that I realized based on your problem:

• The nan is not obtained from a 1/x division, but instead, it is obtained because the density solver failed to find a solution. By default, the density solver uses the previous solution to start the iteration; however, this mixture presents two local minima, and the solver jumps between them; hence, at this point, the solution of the previous solution is very off and makes the solver fail. I added a few extra lines that check if there is a nan either in the density solver solution and if that's the case, the calculation is repeated without relying on any initial guess (I have a routine to automatically generate an initial guess, it is usually more expensive to run than using the previous solution as an initial guess).
• The second problem arises from the amine concentration in the vapor phase. The solution finds that the molar fraction of AEEA is of the order of 1e-10. This extremely small concentration cause big numerical issues because is of the order of the tolerances. One "easy" solution would be to prevent the amine to be in the vapor phase. The predicted boiling point of AEEA is ~500K so that assumptions wouldn't be so crazy.

See some solutions with the above (I'm preventing the ASS step because is jumping between solutions)

• Considering AEEA in the vapor phase:

sol=flash(x0, y0, 'LV', z, T, P, eos, minimization_approach='helmholtz', 
          not_in_x_list=[], not_in_y_list=[], full_output=True, nacc=0) # Set the amount of N2 in the liquid phase to zero
      T: 313.15
      P: 101000
   beta: 0.8175545263742379
  error: 0.025054419352129573
   iter: 75
      X: array([1.47176656e-04, 3.79806416e-02, 6.86335850e-02, 8.93238597e-01])
     v1: 2.4101360142835413e-05
  Xass1: array([4.47088888e-02, 9.72595681e-03, 1.74872530e-01, 6.97966918e-03,
       3.66694469e-01, 2.86091422e-05, 6.01433402e-01, 4.42699813e-02,
       8.89444820e-02, 9.45413699e-02])
 state1: 'L'
      Y: array([8.00558638e-01, 1.38049560e-01, 3.91261490e-10, 6.13918012e-02])
     v2: 0.025721154626608105
  Xass2: array([0.99952635, 1.        , 0.99783144, 0.05341241, 0.99920981,
       0.00156895, 0.99981808, 0.99036564, 0.99641074, 0.9958782 ])
 state2: 'V'
 method: 'Helmholtz_minimization'

• Without considering AEEA in the vapor phase:

sol=flash(x0, y0, 'LV', z, T, P, eos, minimization_approach='helmholtz', 
          not_in_x_list=[], not_in_y_list=[2], full_output=True, nacc=0) # Set the amount of N2 in the liquid phase to zero

      T: 313.15
      P: 101000
   beta: 0.8175545274048039
  error: 2.800413437049716e-08
   iter: 27
      X: array([1.47176607e-04, 3.79806427e-02, 6.86335872e-02, 8.93238594e-01])
     v1: 2.4101360329613608e-05
  Xass1: array([4.47088873e-02, 9.72595645e-03, 1.74872531e-01, 6.97966927e-03,
       3.66694470e-01, 2.86091425e-05, 6.01433400e-01, 4.42699816e-02,
       8.89444821e-02, 9.45413704e-02])
 state1: 'L'
      Y: array([0.80055864, 0.13804956, 0.        , 0.0613918 ])
     v2: 0.02572115464609747
  Xass2: array([0.99952635, 1.        , 0.99783144, 0.05341241, 0.99920981,
       0.00156895, 0.99981808, 0.99036564, 0.99641074, 0.9958782 ])
 state2: 'V'
 method: 'Helmholtz_minimization'

I'll update a new version that checks if the density solver converged correctly shortly.

Edit: I just updated a new version to pypi. For now, I only uploaded a precompiled version for OSX; soon, I'll upload a compiled version for Windows.

Regards,
Gustavo

[YouHyin]

First of all, I'm really sorry for the late response. Thank you so much for your response and updating the package so quickly. I'll refer to your answer and proceed with the flash calculation. Thanks to you, it was very helpful.
Thank you always.

Best whises!
YouHyin