Benchmark: aqueous LLE (water + n-butyl acetate, water + 1-butanol) with esper2023 + rehner2023_binary, and two tp_flash convergence failures #407

For your information, a follow-up to #406 (not a bug report): the same evaluation for aqueous liquid–liquid equilibria with the recommended esper2023.json + rehner2023_binary.json. The parameters README already says that aqueous systems may need more specialised parametrisations, so the accuracy part is probably expected. We post it as a data point, together with two tp_flash convergence failures that may be of more interest.

Summary

• The ester solubility in water is about 80 times too high.
• For water + 1-butanol, the two-phase region closes at about 300–310 K. The measured upper critical solution temperature is 397.55–398.5 K (IUPAC-NIST SDS vol. 15, critical evaluation).
• Convergence: PhaseEquilibrium.tp_flash raises RuntimeError: `stability analysis` did not converge in two cases (reproducer below):
  • the ternary water + 1-propanol + n-butyl acetate at 323.15 K, with the feed at the midpoint of a measured tie line (Costa et al. 2009);
  • water + 1-butanol at 313.15 K.
    For the ternary, an independent tangent-plane test (Michelsen successive substitution) on State.ln_symmetric_activity_coefficient() of the same equation of state finds the feed unstable, a negative tangent-plane distance, so the model does split there. Our own split iteration also converges slowly for that feed.

Reproducer (run from the root of a v0.10.5 checkout)

"""Aqueous LLE with esper2023 + rehner2023_binary (feos 0.10.5) against measured data.

Run from the root of a feos checkout (v0.10.5): python repro_lle.py
"""
import numpy as np
import feos
import si_units as si

P = 101.325 * si.KILO * si.PASCAL


def eos_of(names):
    return feos.EquationOfState.pcsaft(
        feos.Parameters.from_json(
            names, "parameters/pcsaft/esper2023.json", "parameters/pcsaft/rehner2023_binary.json"
        )
    )


def lle_diagram(names, feed, t_min, t_max, n):
    diagram = feos.PhaseDiagram.lle(
        eos_of(names), P, feed=np.array(feed) * si.MOL,
        min_tp=t_min * si.KELVIN, max_tp=t_max * si.KELVIN, npoints=n,
    )
    for a, b in zip(diagram.liquid, diagram.vapor):
        xa, xb = np.asarray(a.molefracs), np.asarray(b.molefracs)
        aq, org = (xa, xb) if xa[0] > xb[0] else (xb, xa)
        print(f"  {a.temperature / si.KELVIN:7.2f} K  water-rich x_2 {aq[1]:.4f}   organic x_2 {org[1]:.4f}")


print("water (1) + butyl acetate (2), PhaseDiagram.lle at 1 atm")
lle_diagram(["water", "butyl acetate"], [0.6, 0.4], 293.15, 323.15, 4)
print("water (1) + 1-butanol (2), PhaseDiagram.lle at 1 atm")
lle_diagram(["water", "1-butanol"], [0.75, 0.25], 288.15, 318.15, 4)

print("tp_flash at 1 atm (feed between the measured liquids)")
for names, t, feed in (
    (["water", "butyl acetate"], 323.15, [0.565, 0.435]),
    (["water", "1-propanol", "butyl acetate"], 323.15, [0.592, 0.091, 0.317]),
    (["water", "1-butanol"], 313.15, [0.75, 0.25]),
):
    try:
        eq = feos.PhaseEquilibrium.tp_flash(eos_of(names), t * si.KELVIN, P, np.array(feed))
        print(f"  {names} {t} K: {np.round(eq.liquid.molefracs, 4)} / {np.round(eq.vapor.molefracs, 4)}")
    except RuntimeError as failure:
        print(f"  {names} {t} K: RuntimeError: {failure}")

Output (feos 0.10.5, si-units 0.11.5, Python 3.15, macOS arm64)

water (1) + butyl acetate (2), PhaseDiagram.lle at 1 atm
   293.15 K  water-rich x_2 0.0055   organic x_2 0.9299
   303.15 K  water-rich x_2 0.0080   organic x_2 0.8981
   313.15 K  water-rich x_2 0.0116   organic x_2 0.8559
   323.15 K  water-rich x_2 0.0166   organic x_2 0.8020
water (1) + 1-butanol (2), PhaseDiagram.lle at 1 atm
   288.15 K  water-rich x_2 0.0304   organic x_2 0.5618
   298.15 K  water-rich x_2 0.0497   organic x_2 0.4501
tp_flash at 1 atm (feed between the measured liquids)
  ['water', 'butyl acetate'] 323.15 K: [0.9834 0.0166] / [0.198 0.802]
  ['water', '1-propanol', 'butyl acetate'] 323.15 K: RuntimeError: `stability analysis` did not converge within the maximum number of iterations.
  ['water', '1-butanol'] 313.15 K: RuntimeError: `stability analysis` did not converge within the maximum number of iterations.

Data sources (read from the NIST TRC ThermoML files unless noted; mass fractions converted with M = 18.015 and 74.123 g/mol)

• H. F. Costa et al., J. Chem. Eng. Data 54 (2009) 2845, doi:10.1021/je900337a, ThermoML data set 27 (binary edge of the ternary data).
• J.-T. Chen, Y.-M. Lin, Fluid Phase Equilib. 259 (2007) 189, doi:10.1016/j.fluid.2007.07.011, data sets 2 and 3.
• Y.-Y. Wu et al., Fluid Phase Equilib. 265 (2008) 1, doi:10.1016/j.fluid.2007.12.003, data set 8 (mass fractions of water 0.92 and 0.203 at 298.15 K).
• V. Gomis et al., J. Chem. Thermodyn. 79 (2014) 242, doi:10.1016/j.jct.2014.07.021, data sets 3 and 4.
• A. F. M. Barton (ed.), IUPAC-NIST Solubility Data Series vol. 15, Alcohols with Water (1984), critical evaluation of 1-butanol + water (UCST 397.55, 397.85, 398.30 and 398.5 K as reported), https://srdata.nist.gov/solubility/IUPAC/SDS-15/SDS-15.pdf.

For context: the published LLE-fitted parametrisation of Veith et al. (ACS Omega 2 (2017) 8982, Fuchs et al. water with σ(T), Gross & Sadowski 1-butanol, k_ij(T) from their Table 2), run in feos by rebuilding the parameters at each temperature, gives 0.013–0.015 / 0.42–0.45 at 288–318 K and 0.0197 / 0.3625 at 365.85 K. That is close to the data at 366 K, but it also does not reproduce the UCST.

🤖 Generated with Claude Code