PhaseEquilibrium.heteroazeotrope_t fails to converge in v0.9 for a case that converges successfully in v0.8 #319
Description
The system is ethane + nitrogen, modelled with PCP-SAFT. Parameters below.
Ethane:
"molarweight": 30.047,
"m": 1.6069,
"sigma": 3.5206,
"epsilon_k": 191.420
Nitrogen (from "gross2005_fit.json"):
"molarweight": 28.01,
"m": 1.1879,
"sigma": 3.3353,
"epsilon_k": 90.99,
"q": 1.1151
BIP: "k_ij": 0.03851
For reference, the three-phase state is:
- xIN2 = 0.3334
- xIIN2 = 0.9683
- yN2 = 0.9984
- P = 3.022 MPa
This is technically not a "heteroazeotrope", just a type III binary with VLLE (vapor is richer in N2 than both liquid phases). However, FeOs's PhaseEquilibrium.heteroazeotrope can still find the three-phase state, and PhaseDiagram.binary_vlle still works, though some manual slicing of the results is necessary. Below is the phase diagram for this system @ T = 125 K, computed with FeOs v0.8:
After upgrading to v0.9, PhaseEquilibrium.heteroazeotrope (specifically, heteroazeotrope_t, when temperature is specified), fails.
My code, v0.8:
from si_units import *
from feos.pcsaft import *
from feos.dft import *
import numpy as np
components = ["ethane", "nitrogen"]
params = PcSaftParameters.from_multiple_json([([components[0]], "..\\nikolaidis2024.json"), ([components[1]], "..\\gross2005_fit.json")], "..\\nikolaidis2024_binary.json")
pcsaft = HelmholtzEnergyFunctional.pcsaft(params)
print(PhaseEquilibrium.heteroazeotrope(pcsaft, 125*KELVIN, (0.03, 0.67), tp_init=3*MEGA*PASCAL, verbosity=Verbosity.Iter, verbosity_bd=Verbosity.Iter))
v0.8 converges successfully:
res outer loop | res inner loop | temperature | molefracs second phase
-------------------------------------------------------------------------------------
| | 3.00000000 MPa | [0.00028902, 0.99971098]
| 3.18940094e-3 | 3.02688871 MPa | [0.00028902, 0.99971098]
| 5.39775033e-5 | 3.02736202 MPa | [0.00028902, 0.99971098]
| 1.72360410e-8 | 3.02736217 MPa | [0.00028902, 0.99971098]
| 2.26485497e-14 | 3.02736217 MPa | [0.00028902, 0.99971098]
4.43572928e0 | | |
| 1.36207899e-5 | 3.02748180 MPa | [0.00157064, 0.99842936]
| 1.10877796e-9 | 3.02748181 MPa | [0.00157064, 0.99842936]
| 1.55431223e-15 | 3.02748181 MPa | [0.00157064, 0.99842936]
2.13205699e-2 | | |
| 9.68720681e-9 | 3.02748189 MPa | [0.00160407, 0.99839593]
| 1.55431223e-15 | 3.02748189 MPa | [0.00160407, 0.99839593]
5.64652750e-4 | | |
1.90669356e-3 | | 3.02748189 MPa | [0.00160500, 0.99839500] NEWTON
| 1.62536651e-13 | 3.02748189 MPa | [0.00160500, 0.99839500]
1.28786648e-10 | | |
4.37585557e-10 | | 3.02748189 MPa | [0.00160500, 0.99839500] NEWTON
2.55801747e-13 | | 3.02748189 MPa | [0.00160500, 0.99839500] NEWTON
Bubble/dew point: calculation converged in 7 step(s)
res outer loop | res inner loop | temperature | molefracs second phase
-------------------------------------------------------------------------------------
| | 3.00000000 MPa | [0.00029465, 0.99970535]
| 2.60481862e-4 | 2.99814232 MPa | [0.00029465, 0.99970535]
| 2.67062532e-7 | 2.99814422 MPa | [0.00029465, 0.99970535]
| 6.23057161e-13 | 2.99814422 MPa | [0.00029465, 0.99970535]
4.24236887e0 | | |
| 1.15998951e-5 | 2.99822680 MPa | [0.00154429, 0.99845571]
| 5.29189803e-10 | 2.99822680 MPa | [0.00154429, 0.99845571]
1.86016636e-2 | | |
| 6.33935593e-9 | 2.99822685 MPa | [0.00157297, 0.99842703]
| 4.22217816e-13 | 2.99822685 MPa | [0.00157297, 0.99842703]
4.32062064e-4 | | |
1.27604271e-3 | | 2.99822685 MPa | [0.00157367, 0.99842633] NEWTON
| 2.07611706e-14 | 2.99822685 MPa | [0.00157367, 0.99842633]
5.65121283e-11 | | |
Bubble/dew point: calculation converged in 5 step(s)
phase 0: T = 125.00000 K, ρ = 5.58275 kmol/m³, x = [0.00163, 0.99837]
phase 1: T = 125.00000 K, ρ = 18.94391 kmol/m³, x = [0.03170, 0.96830]
phase 2: T = 125.00000 K, ρ = 22.35085 kmol/m³, x = [0.66665, 0.33335]
(unrelated, but it appears that the plain verbosity argument (i.e., not verbosity_bd) does nothing here?)
For clarity, due the Python module changes, here is all the code again for v0.9 (using updated parameter file format):
from si_units import *
from feos import *
import numpy as np
components = ["ethane", "nitrogen"]
params = Parameters.from_multiple_json([([components[0]], "..\\nikolaidis2024.json"), ([components[1]], "..\\gross2005_fit.json")], "..\\nikolaidis2024_binary.json")
pcsaft = HelmholtzEnergyFunctional.pcsaft(params)
print(PhaseEquilibrium.heteroazeotrope(pcsaft, 125*KELVIN, (0.03, 0.67), tp_init=3*MEGA*PASCAL, verbosity=Verbosity.Iter, verbosity_bd=Verbosity.Iter))
On v0.9, the same call to PhaseEquilibrium.heteroazeotrope fails:
res outer loop | res inner loop | temperature | pressure | molefracs second phase
--------------------------------------------------------------------------------------------------------
| | 125.00000000 K | 3.00000000 MPa | [0.00020163, 0.69743843]
| 5.39212155e-2 | 125.00000000 K | 2.71451225 MPa | [0.00028902, 0.99971098]
| 4.47494299e-2 | 125.00000000 K | 2.98678125 MPa | [0.00028902, 0.99971098]
| 4.78070819e-3 | 125.00000000 K | 3.02635074 MPa | [0.00028902, 0.99971098]
| 1.15361666e-4 | 125.00000000 K | 3.02736148 MPa | [0.00028902, 0.99971098]
| 7.85498881e-8 | 125.00000000 K | 3.02736217 MPa | [0.00028902, 0.99971098]
4.43572928e0 | | | |
| 1.36207899e-5 | 125.00000000 K | 3.02748180 MPa | [0.00157064, 0.99842936]
| 1.10883991e-9 | 125.00000000 K | 3.02748181 MPa | [0.00157064, 0.99842936]
| 2.36477504e-14 | 125.00000000 K | 3.02748181 MPa | [0.00157064, 0.99842936]
2.13205700e-2 | | | |
| 9.68721392e-9 | 125.00000000 K | 3.02748189 MPa | [0.00160407, 0.99839593]
| 4.10782519e-15 | 125.00000000 K | 3.02748189 MPa | [0.00160407, 0.99839593]
5.64652750e-4 | | | |
| 1.90669356e-3 | 125.00000000 K | 3.02748189 MPa | [0.00160500, 0.99839500] NEWTON
| 1.54458668e-11 | 125.00000000 K | 3.02748189 MPa | [0.00160498, 0.99839502]
1.52780561e-5 | | | |
| 5.16224340e-5 | 125.00000000 K | 3.02748189 MPa | [0.00160500, 0.99839500] NEWTON
| 1.19387760e-11 | 125.00000000 K | 3.02748189 MPa | [0.00160500, 0.99839500] NEWTON
Bubble/dew point: calculation converged in 7 step(s)
res outer loop | res inner loop | temperature | pressure | molefracs second phase
--------------------------------------------------------------------------------------------------------
| | 125.00000000 K | 3.00000000 MPa | [0.00020485, 0.69501935]
| 5.78361502e-2 | 125.00000000 K | 2.71451225 MPa | [0.00029465, 0.99970535]
| 4.79667061e-2 | 125.00000000 K | 2.96717841 MPa | [0.00029465, 0.99970535]
| 4.44839948e-3 | 125.00000000 K | 2.99764634 MPa | [0.00029465, 0.99970535]
| 7.00008855e-5 | 125.00000000 K | 2.99814408 MPa | [0.00029465, 0.99970535]
| 1.91331808e-8 | 125.00000000 K | 2.99814422 MPa | [0.00029465, 0.99970535]
4.24236887e0 | | | |
| 1.15998952e-5 | 125.00000000 K | 2.99822680 MPa | [0.00154429, 0.99845571]
| 5.29242428e-10 | 125.00000000 K | 2.99822680 MPa | [0.00154429, 0.99845571]
1.86016636e-2 | | | |
| 6.33919628e-9 | 125.00000000 K | 2.99822685 MPa | [0.00157297, 0.99842703]
| 7.32747196e-15 | 125.00000000 K | 2.99822685 MPa | [0.00157297, 0.99842703]
4.32062067e-4 | | | |
| 1.27604268e-3 | 125.00000000 K | 2.99822685 MPa | [0.00157367, 0.99842633] NEWTON
| 3.25117711e-12 | 125.00000000 K | 2.99822685 MPa | [0.00157365, 0.99842635]
1.02245376e-5 | | | |
| 3.02099805e-5 | 125.00000000 K | 2.99822685 MPa | [0.00157367, 0.99842633] NEWTON
| 7.33889374e-12 | 125.00000000 K | 2.99822685 MPa | [0.00157367, 0.99842633] NEWTON
Bubble/dew point: calculation converged in 7 step(s)
Traceback (most recent call last):
File "e:\Work\Phase Equilibrium\Phase Diagrams\Mine\binary_vlle_plot.py", line 80, in <module>
print(PhaseEquilibrium.heteroazeotrope(pcsaft, T, (0.03, 0.67), tp_init=3*MEGA*PASCAL, verbosity=Verbosity.Iter, verbosity_bd=Verbosity.Iter))
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
RuntimeError: `PhaseEquilibrium::heteroazeotrope_t` encountered illegal values during the iteration.
Looks like v0.9's first iteration has mole fractions that don't sum to 1? But otherwise both outer loops seem to arrive at the same point. So perhaps something has broken in the state validation afterwards? Some initializations did report invalid (negative) densities.
I tested with various initializations, both closer and farther from the actual values, and could not find an init that converged to the three-phase solution for v0.9; some did converge to a degenerate solution (2/3 phases the same), without reporting a trivial solution as an exception. Also tweaked tolerances and max iter counts up and down, and nothing helped.
I also tested these by specifying the pressure (using the converged value from the v0.8 run, 3.0224492541277184 MPa) rather than temperature, and this time both v0.8 and v0.9 converged successfully (to 124.99999999997081 K)! So the issue is confined to heteroazeotrope_t.