# 16. Phase equilibrium#

It is not necessary to use a Stream object to use phase equilibrium methods. In fact, thermosteam makes it just as easy to compute vapor-liquid equilibrium, bubble and dew points, and fugacities.

## 16.1. Fugacities#

The easiest way to calculate fugacities is through LiquidFugacities and GasFugacities objects:

[1]:

import thermosteam as tmo
import numpy as np
chemicals = tmo.Chemicals(['Water', 'Ethanol'])
tmo.settings.set_thermo(chemicals)

# Create a LiquidFugacities object
F_l = tmo.equilibrium.LiquidFugacities(chemicals)

# Compute liquid fugacities
liquid_molar_composition = np.array([0.72, 0.28])
f_l = F_l(x=liquid_molar_composition, T=355)
f_l

[1]:

array([43338.226, 57731.001])

[2]:

# Create a GasFugacities object
F_g = tmo.equilibrium.GasFugacities(chemicals)

# Compute gas fugacities
gas_molar_composition = np.array([0.43, 0.57])
f_g = F_g(y=gas_molar_composition, T=355, P=101325)
f_g

[2]:

array([43569.75, 57755.25])


## 16.2. Bubble and dew points#

Similarly bubble and dew point can be calculated through BubblePoint and DewPoint objects:

[3]:

# Create a BubblePoint object
BP = tmo.equilibrium.BubblePoint(chemicals)
molar_composition = np.array([0.5, 0.5])

# Solve bubble point at constant temperature
bp = BP(z=molar_composition, T=355)
bp

[3]:

BubblePointValues(T=355.00, P=109407, IDs=('Water', 'Ethanol'), z=[0.5 0.5], y=[0.344 0.656])

[4]:

# Note that the result is a BubblePointValues object which contain all results as attibutes
(bp.T, bp.P, bp.IDs, bp.z, bp.y)

[4]:

(355,
109406.50434728098,
('Water', 'Ethanol'),
array([0.5, 0.5]),
array([0.344, 0.656]))

[5]:

# Solve bubble point at constant pressure
BP(z=molar_composition, P=2*101325)

[5]:

BubblePointValues(T=371.86, P=202650, IDs=('Water', 'Ethanol'), z=[0.5 0.5], y=[0.351 0.649])

[6]:

# Create a DewPoint object
DP = tmo.equilibrium.DewPoint(chemicals)

# Solve for dew point at constant temperautre
dp = DP(z=molar_composition, T=355)
dp

[6]:

DewPointValues(T=355.00, P=92008, IDs=('Water', 'Ethanol'), z=[0.5 0.5], x=[0.849 0.151])

[7]:

# Note that the result is a DewPointValues object which contain all results as attibutes
(dp.T, dp.P, dp.IDs, dp.z, dp.x)

[7]:

(355,
92008.17334352719,
('Water', 'Ethanol'),
array([0.5, 0.5]),
array([0.849, 0.151]))

[8]:

# Solve for dew point at constant pressure
DP(z=molar_composition, P=2*101324)

[8]:

DewPointValues(T=376.25, P=202648, IDs=('Water', 'Ethanol'), z=[0.5 0.5], x=[0.83 0.17])


## 16.3. Vapor liquid equilibrium#

Vapor-liquid equilibrium can be calculated through a VLE object:

[9]:

# First create a material indexer for the VLE object to manage material data
imol = tmo.indexer.MaterialIndexer(l=[('Water', 0.5), ('Ethanol', 0.5)],
g=[('Water', 0.5), ('Ethanol', 0.5)])

# Create a VLE object
vle = tmo.equilibrium.VLE(imol)
vle

[9]:

VLE(imol=MaterialIndexer(
g=[('Water', 0.5), ('Ethanol', 0.5)],
l=[('Water', 0.5), ('Ethanol', 0.5)]),
thermal_condition=ThermalCondition(T=298.15, P=101325))


You can call the VLE object by setting 2 degrees of freedom from the following list:

• T - Temperature [K]

• P - Pressure [Pa]

• V - Molar vapor fraction

• H - Enthalpy [kJ/hr]

• S - Entropy [kJ/K/hr]

• y - Binary molar vapor composition

• x - Binary molar liquid composition

Here are some examples:

[10]:

vle(T=355, P=101325)
vle

[10]:

VLE(imol=MaterialIndexer(
g=[('Water', 0.6361), ('Ethanol', 0.8548)],
l=[('Water', 0.3639), ('Ethanol', 0.1452)]),
thermal_condition=ThermalCondition(T=355.00, P=101325))

[11]:

mixture_enthalpy = vle.mixture.xH(imol, *vle.thermal_condition)
vle(H=mixture_enthalpy, P=202650)
vle

[11]:

VLE(imol=MaterialIndexer(
g=[('Water', 0.6081), ('Ethanol', 0.8183)],
l=[('Water', 0.3919), ('Ethanol', 0.1817)]),
thermal_condition=ThermalCondition(T=373.69, P=202650))

[12]:

vle(V=0.5, P=101325)
vle

[12]:

VLE(imol=MaterialIndexer(
g=[('Water', 0.3874), ('Ethanol', 0.6126)],
l=[('Water', 0.6126), ('Ethanol', 0.3874)]),
thermal_condition=ThermalCondition(T=353.94, P=101325))

[13]:

vle(V=0.5, T=360)
vle

[13]:

VLE(imol=MaterialIndexer(
g=[('Water', 0.3899), ('Ethanol', 0.6101)],
l=[('Water', 0.6101), ('Ethanol', 0.3899)]),
thermal_condition=ThermalCondition(T=360.00, P=127822))

[14]:

vle(x=np.array([0.8, 0.2]), P=101325)
vle

[14]:

VLE(imol=MaterialIndexer(
g=[('Water', 0.8434), ('Ethanol', 0.9609)],
l=[('Water', 0.1566), ('Ethanol', 0.03914)]),
thermal_condition=ThermalCondition(T=356.31, P=127822))

[15]:

vle(y=np.array([0.4, 0.6]), T=360)
vle

[15]:

VLE(imol=MaterialIndexer(
g=[('Water', 0.4627), ('Ethanol', 0.6941)],
l=[('Water', 0.5373), ('Ethanol', 0.3059)]),
thermal_condition=ThermalCondition(T=356.31, P=126587))


Note that some compositions are infeasible; so it is not advised to pass x or y unless you know what you’re doing:

[16]:

vle(x=np.array([0.2, 0.8]), P=101325)
vle

---------------------------------------------------------------------------
InfeasibleRegion                          Traceback (most recent call last)
~\AppData\Local\Temp\ipykernel_14104\2096573133.py in <cell line: 1>()
----> 1 vle(x=np.array([0.2, 0.8]), P=101325)
2 vle

~\Code\biosteam\thermosteam\thermosteam\equilibrium\vle.py in __call__(self, T, P, V, H, S, x, y)
389                         thermal_condition.P = P
390             elif x_spec:
--> 391                 self.set_Px(P, np.asarray(x))
392             else: # y_spec
393                 self.set_Py(P, np.asarray(y))

~\Code\biosteam\thermosteam\thermosteam\equilibrium\vle.py in set_Px(self, P, x)
671         assert self._N == 2, 'number of species in equilibrium must be 2 to specify x'
672         self._thermal_condition.T, y = self._bubble_point.solve_Ty(x, P)
--> 673         self._lever_rule(x, y)
674
675     def set_Ty(self, T, y):

~\Code\biosteam\thermosteam\thermosteam\equilibrium\vle.py in _lever_rule(self, x, y)
653         split_frac = (self._z[0]-x[0])/(y[0]-x[0])
654         if not -0.00001 < split_frac < 1.00001:
--> 655             raise InfeasibleRegion('phase composition')
656         if split_frac > 1:
657             split_frac = 1

InfeasibleRegion: phase composition is infeasible


## 16.4. Liquid-liquid equilibrium#

Liquid-liquid equilibrium can be calculated through a LLE object:

[17]:

tmo.settings.set_thermo(['Water', 'Octane', 'Butanol'])
imol = tmo.indexer.MolarFlowIndexer(
l=[('Water', 304), ('Butanol', 30)],
L=[('Octane', 100)])
lle = tmo.equilibrium.LLE(imol)
lle(T=360)
lle

[17]:

LLE(imol=MolarFlowIndexer(
L=[('Water', 13.35), ('Octane', 99.98), ('Butanol', 25.7)],
l=[('Water', 290.6), ('Octane', 0.02062), ('Butanol', 4.3)]),
thermal_condition=ThermalCondition(T=360.00, P=101325))


Pressure is not a significant factor in liquid-liquid equilibrium, so only temperature is needed.

## 16.5. Vapor-liquid-liquid equilibrium#

For now, the only way to perform vapor-liquid-liquid equilibrium calculations is through Stream objects:

[18]:

tmo.settings.set_thermo(['Water', 'Ethanol', 'Octane'])
s = tmo.Stream('s', Water=100, Ethanol=10, Octane=100)
s.vlle(T=362, P=101325)
s.show()

MultiStream: s
phases: ('L', 'g', 'l'), T: 362 K, P: 101325 Pa
flow (kmol/hr): (g) Water    97
Ethanol  9.69
Octane   45.8
(l) Water    2.98
Ethanol  0.311
Octane   54.2