Validating Flash Calculations

Finding the solution to multiphase equilibrium calculations is challenging and the topic of continuing research. Many commercial packages offer users a great deal of confidence in their answers, but can they be trusted? Thermo can be used to validate the results from other software or identify defects in them.

The following example uses a natural gas mixture two pseudocomponents C7-C16 and C17+. The properties of pure components are taken from Thermo. To do a perfect comparison, the critical properties from other software packages should be substituted into Thermo. This is example S3 from Fonseca-Pérez (2021). The kijs are from Harding and Floudas (2000), and the original pseudocomponents are from Nagarajan, Cullick, and Griewank (1991).

Fonseca-Pérez, R. M., A. Bonilla-Petriciolet, J. C. Tapia-Picazo, and J. E. Jaime-Leal. “A Reconsideration on the Resolution of Phase Stability Analysis Using Stochastic Global Optimization Methods: Proposal of a Reliable Set of Benchmark Problems.” Fluid Phase Equilibria 548 (November 15, 2021): 113180. https://doi.org/10.1016/j.fluid.2021.113180.

Harding, S. T., and C. A. Floudas. “Phase Stability with Cubic Equations of State: Global Optimization Approach.” AIChE Journal 46, no. 7 (July 1, 2000): 1422–40. https://doi.org/10.1002/aic.690460715.

Nagarajan, N. R., A. S. Cullick, and A. Griewank. “New Strategy for Phase Equilibrium and Critical Point Calculations by Thermodynamic Energy Analysis. Part I. Stability Analysis and Flash.” Fluid Phase Equilibria 62, no. 3 (January 1, 1991): 191–210. https://doi.org/10.1016/0378-3812(91)80010-S.

from thermo import *
from scipy.constants import atm
pure_constants = ChemicalConstantsPackage.constants_from_IDs(
    ['methane', 'ethane', 'propane', 'n-butane', 'n-pentane', 'n-hexane'])

pseudos = ChemicalConstantsPackage(Tcs=[606.28,825.67], Pcs=[25.42*atm, 14.39*atm],
                                   omegas=[0.4019, 0.7987], MWs=[140.0, 325.0])
constants = pure_constants + pseudos

properties = PropertyCorrelationsPackage(constants=constants)

T = 353
P = 38500e3
zs = [0.7212, 0.09205, 0.04455, 0.03123, 0.01273, 0.01361, 0.07215, 0.01248]

kijs = [[0.0, 0.002, 0.017, 0.015, 0.02, 0.039, 0.05, 0.09],
[0.002, 0.0, 0.0, 0.025, 0.01, 0.056, 0.04, 0.055],
[0.017, 0.0, 0.0, 0.0, 0.0, 0.0, 0.01, 0.01],
[0.015, 0.025, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
[0.02, 0.01, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
[0.039, 0.056, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
[0.05, 0.04, 0.01, 0.0, 0.0, 0.0, 0.0, 0.0],
[0.09, 0.055, 0.01, 0.0, 0.0, 0.0, 0.0, 0.0]]

eos_kwargs = dict(Tcs=constants.Tcs, Pcs=constants.Pcs, omegas=constants.omegas, kijs=kijs)

gas = CEOSGas(PRMIX, eos_kwargs, HeatCapacityGases=properties.HeatCapacityGases, T=T, P=P, zs=zs)
liq = CEOSLiquid(PRMIX, eos_kwargs, HeatCapacityGases=properties.HeatCapacityGases, T=T, P=P, zs=zs)
liq2 = CEOSLiquid(PRMIX, eos_kwargs, HeatCapacityGases=properties.HeatCapacityGases, T=T, P=P, zs=zs)
phase_list = [gas, liq, liq]


flashN = FlashVLN(constants, properties, liquids=[liq, liq2], gas=gas)
# flashN.PT_SS_TOL = 1e-18
res = flashN.flash(T=T, P=P, zs=zs)
print('There are %s phases present' %(res.phase_count))
print('Mass densities of each liquid are %f and %f kg/m^3' %(res.liquid0.rho_mass(), res.liquid0.rho_mass()))

There are 2 phases present
Mass densities of each liquid are 527.867861 and 527.867861 kg/m^3

import numpy as np
max_fugacity_err = np.max(np.abs(1-np.array(res.liquid0.fugacities())/res.liquid1.fugacities()))
print('The maximum relative difference in fugacity is %e.' %(max_fugacity_err,))

The maximum relative difference in fugacity is 2.773985e-07.