René's URL Explorer Experiment

{"@context":"https://schema.org","@type":"DiscussionForumPosting","headline":"PhaseEquilibrium.heteroazeotrope_t fails to converge in v0.9 for a case that converges successfully in v0.8","articleBody":"The system is ethane + nitrogen, modelled with PCP-SAFT. Parameters below.\nEthane:\n```\n\"molarweight\": 30.047,\n\"m\": 1.6069,\n\"sigma\": 3.5206,\n\"epsilon_k\": 191.420\n```\nNitrogen (from \"gross2005_fit.json\"):\n```\n\"molarweight\": 28.01,\n\"m\": 1.1879,\n\"sigma\": 3.3353,\n\"epsilon_k\": 90.99,\n\"q\": 1.1151\n```\nBIP: `\"k_ij\": 0.03851`\n\n\nFor reference, the three-phase state is:\n- x\u003csup\u003eI\u003c/sup\u003e\u003csub\u003eN\u003csub\u003e2\u003c/sub\u003e\u003c/sub\u003e = 0.3334\n- x\u003csup\u003eII\u003c/sup\u003e\u003csub\u003eN\u003csub\u003e2\u003c/sub\u003e\u003c/sub\u003e = 0.9683\n- y\u003csub\u003eN\u003csub\u003e2\u003c/sub\u003e\u003c/sub\u003e = 0.9984\n- P = 3.022 MPa\n\nThis 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:\n\n\u003cimg src=\"https://github.com/user-attachments/assets/0058bfa8-0a33-4888-a600-73d8ba5b69ef\" height=\"285\"/\u003e\n\n\u003cimg src=\"https://github.com/user-attachments/assets/1af36a7a-9188-4736-a767-dcaed5d4050b\" height=\"285\"/\u003e\n\n\u003cbr/\u003e\u003cbr/\u003e\n\nAfter upgrading to v0.9, `PhaseEquilibrium.heteroazeotrope` (specifically, `heteroazeotrope_t`, when temperature is specified), fails.\n\nMy code, v0.8:\n```\nfrom si_units import *\nfrom feos.pcsaft import *\nfrom feos.dft import *\nimport numpy as np\n\ncomponents = [\"ethane\", \"nitrogen\"]\nparams = PcSaftParameters.from_multiple_json([([components[0]], \"..\\\\nikolaidis2024.json\"), ([components[1]], \"..\\\\gross2005_fit.json\")], \"..\\\\nikolaidis2024_binary.json\")\npcsaft = HelmholtzEnergyFunctional.pcsaft(params)\n\nprint(PhaseEquilibrium.heteroazeotrope(pcsaft, 125*KELVIN, (0.03, 0.67), tp_init=3*MEGA*PASCAL, verbosity=Verbosity.Iter, verbosity_bd=Verbosity.Iter))\n```\n\nv0.8 converges successfully:\n```\nres outer loop | res inner loop |   temperature    | molefracs second phase\n-------------------------------------------------------------------------------------\n               |                |   3.00000000 MPa | [0.00028902, 0.99971098]\n               | 3.18940094e-3  |   3.02688871 MPa | [0.00028902, 0.99971098]\n               | 5.39775033e-5  |   3.02736202 MPa | [0.00028902, 0.99971098]\n               | 1.72360410e-8  |   3.02736217 MPa | [0.00028902, 0.99971098]\n               | 2.26485497e-14 |   3.02736217 MPa | [0.00028902, 0.99971098]\n4.43572928e0   |                |                  |\n               | 1.36207899e-5  |   3.02748180 MPa | [0.00157064, 0.99842936]\n               | 1.10877796e-9  |   3.02748181 MPa | [0.00157064, 0.99842936]\n               | 1.55431223e-15 |   3.02748181 MPa | [0.00157064, 0.99842936]\n2.13205699e-2  |                |                  |\n               | 9.68720681e-9  |   3.02748189 MPa | [0.00160407, 0.99839593]\n               | 1.55431223e-15 |   3.02748189 MPa | [0.00160407, 0.99839593]\n5.64652750e-4  |                |                  |\n1.90669356e-3  |                |   3.02748189 MPa | [0.00160500, 0.99839500] NEWTON\n               | 1.62536651e-13 |   3.02748189 MPa | [0.00160500, 0.99839500]\n1.28786648e-10 |                |                  |\n4.37585557e-10 |                |   3.02748189 MPa | [0.00160500, 0.99839500] NEWTON\n2.55801747e-13 |                |   3.02748189 MPa | [0.00160500, 0.99839500] NEWTON\nBubble/dew point: calculation converged in 7 step(s)\n\nres outer loop | res inner loop |   temperature    | molefracs second phase\n-------------------------------------------------------------------------------------\n               |                |   3.00000000 MPa | [0.00029465, 0.99970535]\n               | 2.60481862e-4  |   2.99814232 MPa | [0.00029465, 0.99970535]\n               | 2.67062532e-7  |   2.99814422 MPa | [0.00029465, 0.99970535]\n               | 6.23057161e-13 |   2.99814422 MPa | [0.00029465, 0.99970535]\n4.24236887e0   |                |                  |\n               | 1.15998951e-5  |   2.99822680 MPa | [0.00154429, 0.99845571]\n               | 5.29189803e-10 |   2.99822680 MPa | [0.00154429, 0.99845571]\n1.86016636e-2  |                |                  |\n               | 6.33935593e-9  |   2.99822685 MPa | [0.00157297, 0.99842703]\n               | 4.22217816e-13 |   2.99822685 MPa | [0.00157297, 0.99842703]\n4.32062064e-4  |                |                  |\n1.27604271e-3  |                |   2.99822685 MPa | [0.00157367, 0.99842633] NEWTON\n               | 2.07611706e-14 |   2.99822685 MPa | [0.00157367, 0.99842633]\n5.65121283e-11 |                |                  |\nBubble/dew point: calculation converged in 5 step(s)\n\nphase 0: T = 125.00000 K, ρ = 5.58275 kmol/m³, x = [0.00163, 0.99837]\nphase 1: T = 125.00000 K, ρ = 18.94391 kmol/m³, x = [0.03170, 0.96830]\nphase 2: T = 125.00000 K, ρ = 22.35085 kmol/m³, x = [0.66665, 0.33335]\n```\n(unrelated, but it appears that the plain `verbosity` argument (i.e., not `verbosity_bd`) does nothing here?)\n\n\u003cbr/\u003e\u003cbr/\u003e\n\nFor clarity, due the Python module changes, here is all the code again for v0.9 (using updated parameter file format):\n```\nfrom si_units import *\nfrom feos import *\nimport numpy as np\n\ncomponents = [\"ethane\", \"nitrogen\"]\nparams = Parameters.from_multiple_json([([components[0]], \"..\\\\nikolaidis2024.json\"), ([components[1]], \"..\\\\gross2005_fit.json\")], \"..\\\\nikolaidis2024_binary.json\")\npcsaft = HelmholtzEnergyFunctional.pcsaft(params)\n\nprint(PhaseEquilibrium.heteroazeotrope(pcsaft, 125*KELVIN, (0.03, 0.67), tp_init=3*MEGA*PASCAL, verbosity=Verbosity.Iter, verbosity_bd=Verbosity.Iter))\n```\n\nOn v0.9, the same call to `PhaseEquilibrium.heteroazeotrope` fails:\n```\nres outer loop | res inner loop |   temperature  |     pressure     | molefracs second phase\n--------------------------------------------------------------------------------------------------------\n               |                | 125.00000000 K |   3.00000000 MPa | [0.00020163, 0.69743843]\n               | 5.39212155e-2  | 125.00000000 K |   2.71451225 MPa | [0.00028902, 0.99971098]\n               | 4.47494299e-2  | 125.00000000 K |   2.98678125 MPa | [0.00028902, 0.99971098]\n               | 4.78070819e-3  | 125.00000000 K |   3.02635074 MPa | [0.00028902, 0.99971098]\n               | 1.15361666e-4  | 125.00000000 K |   3.02736148 MPa | [0.00028902, 0.99971098]\n               | 7.85498881e-8  | 125.00000000 K |   3.02736217 MPa | [0.00028902, 0.99971098]\n4.43572928e0   |                |                |                  |\n               | 1.36207899e-5  | 125.00000000 K |   3.02748180 MPa | [0.00157064, 0.99842936]\n               | 1.10883991e-9  | 125.00000000 K |   3.02748181 MPa | [0.00157064, 0.99842936]\n               | 2.36477504e-14 | 125.00000000 K |   3.02748181 MPa | [0.00157064, 0.99842936]\n2.13205700e-2  |                |                |                  |\n               | 9.68721392e-9  | 125.00000000 K |   3.02748189 MPa | [0.00160407, 0.99839593]\n               | 4.10782519e-15 | 125.00000000 K |   3.02748189 MPa | [0.00160407, 0.99839593]\n5.64652750e-4  |                |                |                  |\n               | 1.90669356e-3  | 125.00000000 K |   3.02748189 MPa | [0.00160500, 0.99839500] NEWTON      \n               | 1.54458668e-11 | 125.00000000 K |   3.02748189 MPa | [0.00160498, 0.99839502]\n1.52780561e-5  |                |                |                  |\n               | 5.16224340e-5  | 125.00000000 K |   3.02748189 MPa | [0.00160500, 0.99839500] NEWTON      \n               | 1.19387760e-11 | 125.00000000 K |   3.02748189 MPa | [0.00160500, 0.99839500] NEWTON      \nBubble/dew point: calculation converged in 7 step(s)\n\nres outer loop | res inner loop |   temperature  |     pressure     | molefracs second phase\n--------------------------------------------------------------------------------------------------------   \n               |                | 125.00000000 K |   3.00000000 MPa | [0.00020485, 0.69501935]\n               | 5.78361502e-2  | 125.00000000 K |   2.71451225 MPa | [0.00029465, 0.99970535]\n               | 4.79667061e-2  | 125.00000000 K |   2.96717841 MPa | [0.00029465, 0.99970535]\n               | 4.44839948e-3  | 125.00000000 K |   2.99764634 MPa | [0.00029465, 0.99970535]\n               | 7.00008855e-5  | 125.00000000 K |   2.99814408 MPa | [0.00029465, 0.99970535]\n               | 1.91331808e-8  | 125.00000000 K |   2.99814422 MPa | [0.00029465, 0.99970535]\n4.24236887e0   |                |                |                  |\n               | 1.15998952e-5  | 125.00000000 K |   2.99822680 MPa | [0.00154429, 0.99845571]\n               | 5.29242428e-10 | 125.00000000 K |   2.99822680 MPa | [0.00154429, 0.99845571]\n1.86016636e-2  |                |                |                  |\n               | 6.33919628e-9  | 125.00000000 K |   2.99822685 MPa | [0.00157297, 0.99842703]\n               | 7.32747196e-15 | 125.00000000 K |   2.99822685 MPa | [0.00157297, 0.99842703]\n4.32062067e-4  |                |                |                  |\n               | 1.27604268e-3  | 125.00000000 K |   2.99822685 MPa | [0.00157367, 0.99842633] NEWTON      \n               | 3.25117711e-12 | 125.00000000 K |   2.99822685 MPa | [0.00157365, 0.99842635]\n1.02245376e-5  |                |                |                  |\n               | 3.02099805e-5  | 125.00000000 K |   2.99822685 MPa | [0.00157367, 0.99842633] NEWTON      \n               | 7.33889374e-12 | 125.00000000 K |   2.99822685 MPa | [0.00157367, 0.99842633] NEWTON      \nBubble/dew point: calculation converged in 7 step(s)\n\nTraceback (most recent call last):\n  File \"e:\\Work\\Phase Equilibrium\\Phase Diagrams\\Mine\\binary_vlle_plot.py\", line 80, in \u003cmodule\u003e        \n    print(PhaseEquilibrium.heteroazeotrope(pcsaft, T, (0.03, 0.67), tp_init=3*MEGA*PASCAL, verbosity=Verbosity.Iter, verbosity_bd=Verbosity.Iter))\n          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\nRuntimeError: `PhaseEquilibrium::heteroazeotrope_t` encountered illegal values during the iteration.\n```\nLooks 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.\n\nI 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.\n\nI 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`.","author":{"url":"https://github.com/ImagineBaggins","@type":"Person","name":"ImagineBaggins"},"datePublished":"2025-11-18T02:03:07.000Z","interactionStatistic":{"@type":"InteractionCounter","interactionType":"https://schema.org/CommentAction","userInteractionCount":4},"url":"https://github.com/319/feos/issues/319"}