CalebBell / thermo

Thermodynamics and Phase Equilibrium component of Chemical Engineering Design Library (ChEDL)
MIT License
594 stars 114 forks source link

Water properties at low temperatures #15

Closed kburns closed 5 years ago

kburns commented 5 years ago

Thank you very much for this package! I'm finding it very useful, but I'm seeing some unexpected results in the properties of water at low temperatures. In particular, fresh water at 1 atm is known to have a density maximum around 4C, with its density decreasing at both lower and higher temperatures. However, if I plot something like [thermo.Chemical('water', T=Ti).rhol for Ti in np.linspace(273.15, 373.15, 1000)], the returned density seems to be monotonically decreasing at all temperatures. Any thoughts?

CalebBell commented 5 years ago

Hi Dr. Keaton, Are you missing the library 'coolprop'? thermo uses it for high-accuracy property calculations of individual chemicals. The bulk of thermo's property calculations are more focused on chemicals with little data available and complex mixture. I've only been free to do so because of its excellent accuracy for ~100 pure compounds.

There's one more issue at play even if you have CoolProp installed. thermo uses an awful lot of caching to speedup its property calculations. One cached method was to select a property calculation method based on the currently specified T, and P; and then for the next calculation, don't go back and see if the method is still the single best method - just keep using it then on until it is not valid. It pos possible to disable caching by setting "thermo.chemical.caching" to False.

Because you picked a starting temperature of 273.15 K, which is 0.02 K below CoolProp's lowest temperature, that method was found to be not valid. If you pick 273.17 K as your starting point, you will see the correct plot, and the density maximum at 4 degrees.


# thermo.chemical.caching = False # Set this to False or use a higher T min
Ts = np.linspace(273.17, 373.15, 1000)
rhos = [thermo.Chemical('water', T=Ti).rhol for Ti in Ts]
plt.plot(Ts, rhos)
plt.show()
print(Ts[rhos.index(max(rhos))]-273.15)

water_density 4.023203203203252

Sincerely, Caleb

kburns commented 5 years ago

Great! Installing coolprop, turning off the cache, and starting at higher temps worked perfectly. Thanks!