Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[SAFT-Gamma-Mie] Question about VLLE for binary mixture #13

Open
YouHyin opened this issue Mar 28, 2024 · 2 comments
Open

[SAFT-Gamma-Mie] Question about VLLE for binary mixture #13

YouHyin opened this issue Mar 28, 2024 · 2 comments

Comments

@YouHyin
Copy link

YouHyin commented Mar 28, 2024

Hello, I'm trying to reproduce the following paper's a binary mixture of water and triethylamine (TEA) VLLE results (Figure.3) using sgtpy.
Using the example code, I have created the following code.

Paper : Lee, Y. S., Galindo, A., Jackson, G., & Adjiman, C. S. (2023). Enabling the direct solution of challenging computer-aided molecular and process design problems: Chemical absorption of carbon dioxide. Computers & Chemical Engineering, 174, 108204.

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

TEA = component(GC={'CH3': 3, 'CH2': 3,'N': 1})
water = component(GC={'H2O':1})

mix = TEA + water

mix.saftgammamie()
eos = saftgammamie(mix)

# phase equilibria conditions
T = 365.  # K
P = 0.101e6  # Pa
z = np.array([0.15, 0.85])

from sgtpy.equilibrium import tpd_min
# initial guess for the composition of the trial phase
w0 = np.array([0.9, 0.1])
x0, tpx0 = tpd_min(w0, z, T, P, eos, stateW = 'L', stateZ = 'L')
# initial guess for the composition of the trial phase
w0, tpw0 = tpd_min(w0, z, T, P, eos, stateW = 'L', stateZ = 'L')
y0 = np.array([0.2, 0.8])
y0, tpv0 = tpd_min(w0, z, T, P, eos, stateW = 'V', stateZ = 'L')

from sgtpy.equilibrium import vlleb
# initial guesses for aqueous, organic and vapor phase composition obtained from tpd minimization
# calculation at fixed temperature
P0 = 1.01325e5  # Pa
vlleb(x0, w0, y0, P0, T, 'T', eos, full_output=True)

# initial guesses for aqueous, organic and vapor phase composition obtained from tpd minimization
# calculation at fixed pressure
T0 =  363.  # K
sol = vlleb(x0, w0, y0, T0, P, 'P', eos, full_output=True)

from sgtpy.equilibrium import bubbleTy, lle
from sgtpy.equilibrium import vlleb
# three phase equilibria computation
sol = vlleb(x0, w0, y0, T0, P, 'P', eos, full_output=True)
Xhaz,  Whaz, Yhaz, Thaz = sol.X, sol.W, sol.Y, sol.T

n = 100

# VLE zone 1
x1 = np.linspace(0, Xhaz[0], n)
XI = np.array([x1, 1-x1])
YI = np.zeros_like(XI)
TI = np.zeros(n)

vxI = np.zeros(n)
vyI = np.zeros(n)
i = 0
T0 = 373.
sol = bubbleTy(XI[:, i], T0, XI[:, i], P, eos, full_output = True)
YI[:, i], TI[i] = sol.Y, sol.T
vxI[i], vyI[i] = sol.v1, sol.v2
for i in range(1, n):
    sol = bubbleTy(YI[:, i-1], TI[i-1], XI[:, i], P, eos, full_output=True, v0=[vxI[i-1], vyI[i-1]])
    YI[:, i], TI[i] = sol.Y, sol.T
    vxI[i], vyI[i] = sol.v1, sol.v2
    
# VLE zone 2
w1 = np.linspace(Whaz[0], 1, n)
XII = np.array([w1, 1-w1])
YII = np.zeros_like(XII)
TII = np.zeros(n)
vxII = np.zeros(n)
vyII = np.zeros(n)

i = 0
sol = bubbleTy(Yhaz, Thaz, XII[:, i], P, eos, full_output = True)
YII[:, i], TII[i] = sol.Y, sol.T
vxII[i], vyII[i] = sol.v1, sol.v2

for i in range(1, n):
    sol = bubbleTy(YII[:, i-1], TII[i-1], XII[:, i], P, eos, full_output = True, v0=[vxII[i-1], vyII[i-1]])
    YII[:, i], TII[i] = sol.Y, sol.T
    vxII[i], vyII[i] = sol.v1, sol.v2

# LLE calculation
Tll =  np.linspace(Thaz,  290, n)
Xll = np.zeros([2, n])
Wll = np.zeros([2, n])

vxll = np.zeros(n)
vwll = np.zeros(n)

i = 0 
Z = (Xhaz+Whaz)/2
sol = lle(Xhaz, Whaz, Z, Tll[i], P, eos, full_output=True)
Xll[:, i], Wll[:, i] = sol.X
vxll[i], vwll[i] = sol.v

for i in range(1, n):
    Z = (Xll[:, i-1] + Wll[:, i-1])/2
    sol = lle(Xll[:, i-1], Wll[:, i-1], Z, Tll[i], P, eos, full_output=True, v0=[vxll[i-1], vwll[i-1]])
    Xll[:, i], Wll[:, i] = sol.X
    vxll[i], vwll[i] = sol.v

fig = plt.figure(figsize=(5,5))
ax = fig.add_subplot(111)
ax.plot(XI[0], TI, color = 'b')
ax.plot(YI[0], TI, color = 'b')

ax.plot(XII[0], TII, color = 'b')
ax.plot(YII[0], TII, color = 'b')

ax.plot(Xll[0], Tll, color = 'b')
ax.plot(Wll[0], Tll, color = 'b')

ax.plot([Xhaz[0], Yhaz[0], Whaz[0]], [Thaz, Thaz, Thaz], 'o-', color = 'b')


ax.tick_params(direction='in')
ax.set_xlabel('$x_1, y_1$')
ax.set_ylabel('T / K')
ax.set_ylim([280, 400])
ax.set_xlim([0, 1])

Result :
image
The figure I want to reproduce :
image

As a result, there seems to be a problem with the LLE prediction.
What would be the best way to set z and the values of w0, z0, and y0 in this code?
Thanks as always for your help!
I look forward to your response.

@gustavochm
Copy link
Owner

gustavochm commented Apr 2, 2024

Hi YouHyin,

I think there is a problem in your initial guesses for the VLLE. If you get the VLLE correctly you can easily construct the whole diagram. See below.

Edit: Initializing the VLLE is not trivial. Ideally you want to have z what is between the VLLE line. Then you can minimize the tpd to obtain possible initial guesses.


import numpy as np
import matplotlib.pyplot as plt
from sgtpy import component, mixture, saftgammamie
from sgtpy.equilibrium import vlleb, bubbleTy, lle
from sgtpy.equilibrium import tpd_minimas


TEA = component(GC={'CH3': 3, 'CH2': 3,'N': 1})
water = component(GC={'H2O':1})

mix = TEA + water

mix.saftgammamie()
eos = saftgammamie(mix)

P = 0.101e6 # Pa

# obtaining three phase equilibria
T0 = 350. # K

z = np.array([0.6, 0.4])
x0s, tpd0s = tpd_minimas(2, z, T0, P, eos, stateW='L', stateZ='L')
y0s, tpd0 = tpd_minimas(1, z, T0, P, eos, stateW='V', stateZ='L')

x0 = x0s[0]
w0 = x0s[1]
y0 = y0s[0]

sol_vlle = vlleb(x0, w0, y0, T0, P, 'P', eos, full_output=True)
Xhaz,  Whaz, Yhaz, Thaz = sol_vlle.X, sol_vlle.W, sol_vlle.Y, sol_vlle.T

#######

n = 100

# VLE zone 1
x1 = np.linspace(0, Xhaz[0], n)
XI = np.array([x1, 1-x1])
YI = np.zeros_like(XI)
TI = np.zeros(n)

vxI = np.zeros(n)
vyI = np.zeros(n)
i = 0
T0 = 373.
sol = bubbleTy(XI[:, i], T0, XI[:, i], P, eos, full_output = True)
YI[:, i], TI[i] = sol.Y, sol.T
vxI[i], vyI[i] = sol.v1, sol.v2
for i in range(1, n):
    sol = bubbleTy(YI[:, i-1], TI[i-1], XI[:, i], P, eos, full_output=True, v0=[vxI[i-1], vyI[i-1]])
    YI[:, i], TI[i] = sol.Y, sol.T
    vxI[i], vyI[i] = sol.v1, sol.v2

# VLE zone 2
w1 = np.linspace(Whaz[0], 1, n)
XII = np.array([w1, 1-w1])
YII = np.zeros_like(XII)
TII = np.zeros(n)
vxII = np.zeros(n)
vyII = np.zeros(n)

i = 0
sol = bubbleTy(Yhaz, Thaz, XII[:, i], P, eos, full_output = True)
YII[:, i], TII[i] = sol.Y, sol.T
vxII[i], vyII[i] = sol.v1, sol.v2

for i in range(1, n):
    sol = bubbleTy(YII[:, i-1], TII[i-1], XII[:, i], P, eos, full_output = True, v0=[vxII[i-1], vyII[i-1]])
    YII[:, i], TII[i] = sol.Y, sol.T
    vxII[i], vyII[i] = sol.v1, sol.v2

# LLE calculation
Tll =  np.linspace(Thaz,  310, n)
Xll = np.zeros([2, n])
Wll = np.zeros([2, n])

vxll = np.zeros(n)
vwll = np.zeros(n)

i = 0 
Z = (Xhaz+Whaz)/2
sol = lle(Xhaz, Whaz, Z, Tll[i], P, eos, full_output=True)
Xll[:, i], Wll[:, i] = sol.X
vxll[i], vwll[i] = sol.v

for i in range(1, n):
    Z = (Xll[:, i-1] + Wll[:, i-1])/2
    sol = lle(Xll[:, i-1], Wll[:, i-1], Z, Tll[i], P, eos, full_output=True, v0=[vxll[i-1], vwll[i-1]])
    Xll[:, i], Wll[:, i] = sol.X
    vxll[i], vwll[i] = sol.v

Screenshot 2024-04-02 at 13 54 52

I hope that helps,

Gustavo

@YouHyin
Copy link
Author

YouHyin commented Apr 4, 2024

Thank you so much. Thanks to your advice and code, it helped me a lot.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants