Negative vapor and unsteady liquid quality caluculation for Helium

[PerplexedFox]

Hello!

I've noticed some issues during the calculation of Helium state variables. For example, the quality of the pure vapor state results in -1. For a liquid state, it appears unsteady and oscillates between 0 and -1. In comparison, REFPROP (10.0) shows steady curves of 1→100 (superheated vapor) and 0 (saturated liquid). Here are the data points to compare in MATLAB:

p = [137461.323112938;138230.195483501;139009.601323972;139802.273201606;140612.019520862;140624.831673310;140624.881557153];
rho_l = [116.801408146380;116.622690424081;116.440870994868;116.255265192642;116.064936888442;116.061915920152;116.061904153613];
rho_v = [21.1913767379303;21.2996648399277;21.4066320262526;21.5120163010261;21.6154385418442;21.6170372666172;21.6170434816964 ];
 for i = 1:7
q_v_cool(i) = py.CoolProp.CoolProp.PropsSI('Q','P',p(i),'D',rho_v(i),"HELIUM");
q_v_ref(i) = py.CoolProp.CoolProp.PropsSI('Q','P',p(i),'D',rho_v(i),"REFPROP::HELIUM");
q_l_ref(i) = py.CoolProp.CoolProp.PropsSI('Q','P',p(i),'D',rho_l(i),"REFPROP::HELIUM");
q_l_cool(i) = py.CoolProp.CoolProp.PropsSI('Q','P',p(i),'D',rho_l(i),"HELIUM");
end
figure
hold on
plt_v = plot([1;2;3;4;5;6;7],q_v_cool,"DisplayName","Vapor, Coolprop","LineWidth",1.5,"Color","r");
plt_v_r = plot([1;2;3;4;5;6;7],q_v_ref,"LineStyle",'--',"DisplayName","Vapor, REFPROP","LineWidth",1.5,"Color","m");
plt_l = plot([1;2;3;4;5;6;7],q_l_cool,"DisplayName","Liquid, Coolprop","LineWidth",1.5,"Color","b");
plt_l_r = plot([1;2;3;4;5;6;7],q_l_ref,"LineStyle",'--',"DisplayName","Liquid, REFPROP","LineWidth",1.5,"Color","k");
legend

If you repeat it for the enthalpy, you will see that the enthalpy calculated with CoolProp deviates by almost a factor of 10² J/kg.
Is it a bug? Any help would be appreciated!

[ibell]

Where are your state points? Please provide a runnable example with all necessary imports and so on.

[PerplexedFox]

please see the post above.
Liquid density: rho_l = [116.801408146380;116.622690424081;116.440870994868;116.255265192642;116.064936888442;116.061915920152;116.061904153613];
Vapor density rho_v = [21.1913767379303;21.2996648399277;21.4066320262526;21.5120163010261;21.6154385418442;21.6170372666172;21.6170434816964 ];
Pressure: p = [137461.323112938;138230.195483501;139009.601323972;139802.273201606;140612.019520862;140624.831673310;140624.881557153];

[ibell]

Doesn't answer my question: where in the phase diagram (in words) are the state points?

[PerplexedFox]

Hi, sorry for the late answer. The data are for a saturated liquid and a superheated vapor

[ibell]

Looks like a weird bug in CoolProp. @msaitta-mpr can you take a look? The liquid points are basically saturated, the vapor ones are superheated. Of course you could calculate the vapor quality yourself as a workaround.

import matplotlib.pyplot as plt, numpy as np
import CoolProp.CoolProp as CP 

p = [137461.323112938,138230.195483501,139009.601323972,139802.273201606,140612.019520862,140624.831673310,140624.881557153]
rho_l = [116.801408146380,116.622690424081,116.440870994868,116.255265192642,116.064936888442,116.061915920152,116.061904153613]
rho_v = [21.1913767379303,21.2996648399277,21.4066320262526,21.5120163010261,21.6154385418442,21.6170372666172,21.6170434816964 ]

q_v_cool = CP.PropsSI('Q','P',p,'D',rho_v,"HELIUM")
q_v_ref = CP.PropsSI('Q','P',p,'D',rho_v,"REFPROP::HELIUM")
q_l_ref = CP.PropsSI('Q','P',p,'D',rho_l,"REFPROP::HELIUM")
q_l_cool = CP.PropsSI('Q','P',p,'D',rho_l,"HELIUM")

Tt, Tc = [CP.PropsSI(k,'Helium') for k in ['Ttriple','Tcrit']]
Ts = np.linspace(Tt, Tc, 1000)
ps = CP.PropsSI('P','T',Ts,'Q',0,'REFPROP::Helium')
dLs = CP.PropsSI('D','T',Ts,'Q',0,'REFPROP::Helium')
dVs = CP.PropsSI('D','T',Ts,'Q',1,'REFPROP::Helium')
plt.plot(dLs, ps)
plt.plot(dVs, ps)
plt.plot(rho_l, p, 'x')
plt.plot(rho_v, p, 'o')
plt.gca().set(xlabel=r'$\rho$ / kg/m$^3$', ylabel='$p$ / Pa')
plt.tight_layout(pad=0.2)
plt.show()

plt.plot(p, q_v_cool, label='CoolProp')
plt.plot(p, q_v_ref, label='REFPROP')
plt.title('gas (should be > 1)')
plt.legend(loc='best')
plt.gca().set(xlabel='$p$ / Pa', ylabel='Q')
plt.show()

plt.plot(p, q_l_cool, label='CoolProp')
plt.plot(p, q_l_ref, label='REFPROP')
plt.title('liquid (should be = 0)')
plt.legend(loc='best')
plt.gca().set(xlabel='$p$ / Pa', ylabel='Q')
plt.show()

yields

and

and

I am really not sure what is causing this. Someone will have to dig in.