Saturation curves not converging correctly at low temperatures

[GustavoRodriguezM]

Hi! First of all, thanks a lot for all the work on Clapeyron, it’s been super helpful for my project.

I’ve been using it to generate large datasets of saturation and bulk thermodynamic properties for many equations of state. One issue I keep running into is that at low temperatures the saturated vapor/liquid densities often come back as NaN.

As shown above, a large fraction of the data at low temperatures is NaN. This behavior appears across thousands of model–compound combinations and is not limited to a specific equation-of-state family, it occurs in cubic, SAFT-based, and group-contribution models alike.

In principle, these states should be recoverable by supplying suitable initial guesses. In my case, I fitted Wagner equations to saturated liquid and vapor densities, as well as vapor-pressure data, achieving very low fitting errors. Also, notice that the vapor pressure itself is computed correctly, while the corresponding saturated densities come out wrong. Even after applying these fitted correlations as initial guesses, the solver fails to recover the correct densities: saturated liquid densities collapse toward zero, and vapor densities diverge.

This behavior appears to be linked to incorrect phase identification in parts of the phase diagram. Below the first temperature at which failures occur, part of the liquid phase is sometimes misidentified and assigned densities that fall within the vapor density range.

For reference, the code used to generate all data is provided below.

using Clapeyron, PyCall
scipy = pyimport("scipy")
plt = pyimport("matplotlib.pyplot")
np = pyimport("numpy")
least_squares = pyimport("scipy.optimize").least_squares

function vapor_pressure(T, Tc, Pc, param_values)
    tau = 1 .- T ./ Tc
    params = vcat(collect(param_values), zeros(4 - length(param_values)))
    ln_p = (Tc ./ T) .* (params[1] .* tau .+ params[2] .* tau.^1.5 .+ params[3] .* tau.^2 .+ params[4] .* tau.^4)
    return Pc .* np.exp.(ln_p)
end

function rho_liquid(T, Tc, rho_c, param_values)
    tau = 1 .- T ./ Tc
    params = vcat(collect(param_values), zeros(4 - length(param_values)))
    ln_rho = (params[1] .* tau.^(2/6) .+ params[2] .* tau.^(3/6) .+ params[3] .* tau.^(7/6) .+ params[4] .* tau.^(9/6))
    return rho_c .* np.exp.(ln_rho)
end

function rho_vapor(T, Tc, rho_c, param_values)
    tau = 1 .- T ./ Tc
    params = vcat(collect(param_values), zeros(5 - length(param_values)))
    ln_rho = (params[1] .* tau.^(2/6) .+ params[2] .* tau.^(4/6) .+ params[3] .* tau.^(7/6) .+
              params[4] .* tau.^(13/6) .+ params[5] .* tau.^(25/6))
    return rho_c .* np.exp.(ln_rho)
end

function adjuster(x, y, fun, N, method="trf")
    params = []
    AAD = 1e7

    for num_param in 1:N
        guess = ones(num_param)
        if length(params) > 0
            guess[1:length(params)] .= params
        end
        # residuals to minimize
        function residuals(p)
            return fun(x, p...) - y
        end
        res = least_squares(residuals, guess, method=method, max_nfev=10000)
        params = res["x"]
        pred = fun(x, params...)
        AAD = np.mean(np.abs(pred - y)*100/y)
    end
    return params, AAD
end

function fit_densities(Tc, rho_c, T, rho_sat_vapor, rho_sat_liquid)

    f_liq,f_vap,liquid_params,vapor_params = nothing, nothing, nothing, nothing

    wrap_liq(T, p...) = rho_liquid(T, Tc, rho_c, p)
    wrap_vap(T, p...) = rho_vapor(T, Tc, rho_c, p)

    mask_below_Tc = T .< Tc
    mask_T_min = T .> Tc * 0
    is_not_nan = np.isfinite(rho_sat_liquid) .& np.isfinite(rho_sat_vapor)

    T = T[mask_below_Tc][mask_T_min][is_not_nan]
    rho_sat_liquid = rho_sat_liquid[mask_below_Tc][mask_T_min][is_not_nan]
    rho_sat_vapor = rho_sat_vapor[mask_below_Tc][mask_T_min][is_not_nan]

    liquid_params,AAD_liq = adjuster(T,rho_sat_liquid,wrap_liq,4)
    f_liq = T -> wrap_liq(T, liquid_params...)

    vapor_params,AAD_vap = adjuster(T,rho_sat_vapor,wrap_vap,5)
    f_vap = T -> wrap_vap(T, vapor_params...)
    
    return f_liq, f_vap, AAD_liq, AAD_vap
end

function fit_vapor_pressure(Tc, Pc, T, P_sat)

    AAD,params,Psat_fun = nothing, nothing, nothing
    wrapper(T, p...) = vapor_pressure(T, Tc, Pc, p)

    mask_below_Tc = T .< Tc    
    mask_T_min = T .> Tc * 0
    is_not_nan = np.isfinite(P_sat)

    P_sat = P_sat[mask_below_Tc][mask_T_min][is_not_nan]
    T = T[mask_below_Tc][mask_T_min][is_not_nan]

    P_sat = np.append(P_sat, [Pc])
    T = np.append(T, [Tc])
    
    params,AAD = adjuster(T,P_sat,wrapper,4)
    Psat_fun(T) = wrapper(T, params...)

    return Psat_fun,AAD
end


compound = "1-Butene"
CES = "PR" 
N = 500
model = eval(Meta.parse("$(CES)([\"$(compound)\"])"))

(Tc, Pc, Vc) = crit_pure(model)
Tmax = Tc
Tmin = 0.2*Tc
T = LinRange(Tmin, Tmax, N)



Tsat = T
Rho_sat_liq_CES, Sres_sat_liq_CES, Rho_sat_vap_CES, Sres_sat_vap_CES, pv_CES = zeros(length(Tsat)), zeros(length(Tsat)),zeros(length(Tsat)), zeros(length(Tsat)), zeros(length(Tsat))
for (t_idx, t) in enumerate(Tsat)
    (pv, vl, vv) = saturation_pressure(model, t)

    pv_CES[t_idx] = 1 / vl
    Rho_sat_liq_CES[t_idx] = 1 / vl
    Rho_sat_vap_CES[t_idx] = 1 / vv
    Sres_sat_liq_CES[t_idx] = Clapeyron.VT_entropy_res(model, vl, t, [1.])
    Sres_sat_vap_CES[t_idx] = Clapeyron.VT_entropy_res(model, vv, t, [1.])
end




nan_index_liq = Set(np.where(np.isnan(Rho_sat_liq_CES))[1])
nan_index_vap = Set(np.where(np.isnan(Rho_sat_vap_CES))[1])
nan_index_vp = Set(np.where(np.isnan(pv_CES))[1])
nan_index = union(nan_index_liq, nan_index_vap, nan_index_vp)


P_fit, AARD_vp = fit_vapor_pressure(Tc, Pc, Tsat, pv_CES)
liq, vap, AARD_liq, AARD_vap = fit_densities(Tc, 1/Vc, Tsat, Rho_sat_vap_CES, Rho_sat_liq_CES) 

println(AARD_vp, " " ,AARD_liq, "  " , AARD_vap)

for i in nan_index
    ind = i + 1
    t = Tsat[ind]
    
    (pv, vl, vv) = saturation_pressure(model, t, IsoFugacitySaturation(p0 = P_fit(t), vl = liq(t), vv = vap(t), max_iters = 1000, p_tol = 1e-1000 ))

    if isnan(vv) && AARD_vap<=1
        vv = vap(t)
    end

    if isnan(vl) && AARD_liq<=1 
        vl = liq(t)
    end

    if isnan(pv) && AARD_vap<=1
        pv = P_fit(t)
    end
    Rho_sat_liq_CES[ind] = 1 / vl
    Rho_sat_vap_CES[ind] = 1 / vv

    Sres_sat_liq_CES[ind] = Clapeyron.VT_entropy_res(model, vl, t, [1.])
    Sres_sat_vap_CES[ind] = Clapeyron.VT_entropy_res(model, vv, t, [1.])

    pv_CES[ind] = pv
    end

plt.figure()
plt.plot(Tsat ./Tc ,Rho_sat_liq_CES,label="sat liq")
plt.plot(Tsat ./Tc ,Rho_sat_vap_CES,label="sat vap")
plt.xlabel("Tr [-]")
plt.ylabel("Saturation Densities [mol/m^3]")
plt.title("$(CES) for $(compound)")
plt.show()

[simonmb]

We tried using big to have a higher numerical precision on the calls (as recommended to me by @pw0908), but we did not see any improvement. Since the vapor pressure is correct, we also tried using the volume function for each phase with higher numerical precision, but this also did not work.

[longemen3000]

There were some floating point issues on the cubic polynomial solver (that i hope are fixed). for the polynomial volume solver, there is a hard limit on Float64 accuracy, that corresponds to B = b*p*RT⁻¹ > eps(Float64). if B is lower than this threshold, the volume solver will diverge no matter what.

I also added a fix on the initial point function for low temperatures (Clapeyron.x0_sat_pure_near0) that discards a vapour-like liquid solution if it appears. With those two fixes, i can reach Tsat = 0.183Tc on 1-butene on Float64. For BigFloat, setting atol = 1e-12 (over the default atol = 1e-10), allowed me to reach Tsat = 0.05Tc. after that, numerical issues starting appearing with extended precision too.

Note that those changes only take in consideration cubic equations of state. each EoS could have similar floating-point problems.

[simonmb]

@longemen3000, thanks for the fix. @GustavoRodriguezM, correct me if I am wrong, but we also saw this for SAFT based EOSs, right? Would there be a similar fix for these?

[GustavoRodriguezM]

Yes indeed, PCSAFT seems to be the eos with the highest percentage or errors like this, PCPSAFT is likewise even more problematic than the cubics. Other SAFTs show similar behavior

[longemen3000]

For the case of 1-butene error with PCSAFT at T = 0.25Tc

there is an additional spinodal at high pressures:

our problem is that the liquid volume initial point encounters this artifact before reaching the physically valid part of the P-V curve. if i remember correctly, this is a known problem with PCSAFT. But on practical terms, we need a better liquid initial point for PCSAFT (and PCSAFT derivatives)

[longemen3000]

at low temperatures, the mechanically unstable volume at zero pressure can be approximated by $-B_2(model,T,z)$ ,so i just slapped a naive initial point finder for PCSAFT models that kicks in if T < 0.9T_scale(model,z). because the volume solver results cascades through all equilibria procedures, the improved initialization seemingly fixed the saturation pressure issues. I now can get saturation_pressure(model,0.18Tc) in Float64 arithmetic