Wrong values for temperature-dependent sigma of water in PC-SAFT

[pw0908]

The parameters used for the sigma of water here is not the default one from Cameretti and Sadowski (they seem to have been obtained from a paper by Baird et al.):

CoolProp/src/Backends/PCSAFT/PCSAFTFluid.cpp

Lines 58 to 66 in ffa5a27

	void PCSAFTFluid::calc_water_sigma(double t) {
	if (t > 473.16) {
	throw ValueError("The current function for sigma for water is only valid for temperatures below 473.15 K.");
	} else if (t < 273) {
	throw ValueError("The current function for sigma for water is only valid for temperatures above 273.15 K.");
	}
	params.sigma = 3.8395 + 1.2828 * exp(-0.0074944 * t) - 1.3939 * exp(-0.00056029 * t);
	}

The actual values should be:

params.sigma = 2.7927 + 10.1100*exp(-0.01775*t)-1.41700*exp(-0.01146*t)

[ibell]

Zach Baird (@zmeri ), do you have a comment here? Perhaps we need to update the documentation that you implemented your parameters rather than the canonical ones?

[zmeri]

Yes, for our article we needed to calculate properties at higher temperatures, so the original values given by Cameretti and Sadowski did not work because, if I remember correctly, they were only fit for data up to about 100 C. I remember comparing the results of both equations and for lower temperatures they gave the same result to within the accuracy of the PC-SAFT model, if my memory isn't failing me. Are you seeing any significant differences in the results when using ePC-SAFT with our temperature dependent equation for sigma?

[pw0908]

I've actually been doing benchmarking for ePC-SAFT between the pcsaft and CoolProp implementations (for my own implementation in Clapeyron.jl). They don't agree to machine accuracy. The Clapeyron.jl implementation agrees to machine precision with CoolProp though, when using the same temperature-dependent segment diameter.

I also believe the example given in the pcsaft implementation uses the standard correlation?

On that note: how certain are you that your implementations reproduce the literature results? Particularly from the ePC-SAFT revised paper. Despite agreeing to machine precision with CoolProp, I can't reproduce the paper's results and I couldn't find a way to reproduce the paper's results from CoolProp.

[zmeri]

pcsaft and CoolProp use essentially the same code for ePC-SAFT. I wrote the pcsaft code and then copied it over to CoolProp and made the necessary changes to fit it into the CoolProp architecture. I tested the CoolProp implementation against pcsaft, so they should give quite similar results. Are you testing for machine accuracy with values that do not depend on the density? I have before found that differences can come because of differences in when the solvers decide they have achieved an acceptable value and stop. If I remember correctly, there is also an iteration loop to find a value for the association term, so it can also be harder to get exactly equal results (to machine accuracy) for mixtures with the association term.

For both pcsaft and CoolProp I wrote unit tests to check that literature results are reproduced and these are included in the code on Github. Given though that Held and his team have not released their own ePC-SAFT code, I had no way to check ePC-SAFT to machine accuracy, so in the tests I checked they were within expected error limits to the literature values given in the ePC-SAFT articles. I did get some data from Held to check a single datapoint, and I remember that some differences came because he used more decimal places for parameters than were reported in the article. For instance, I remember that in one case a value was given as 3.60 (maybe a sigma value) in the article, but Held used something like 3.5995 in the calculation he sent. Small differences like this don't significantly affect the results, but could cause differences greater than machine accuracy.

[pw0908]

These are the exact benchmarks I ran, along with the conditions, parameters and outputs. I compared a_res, Z and the fugacity coefficients (I couldn't obtain them with CoolProp). The agreement is pretty much spot on between CoolProp and Clapeyron, but there are some discrepancies with pcsaft.
benchmarks_epcsaft.xlsx

Despite getting this agreement with CoolProp and pcsaft, I cannot, for the life of me, reproduce the results from Held's 2014 paper, specifically the osmotic coefficients. Since one can't obtain fugacity coefficients from CoolProp, I couldn't obtain the osmotic coefficients. As for pcsaft, the osmotic coefficient solver seems broken? I did get in touch with Held and I am waiting to hear back regarding benchmarks.

Addressing the original issue here (@ibell): I still think, since it's more widely used, the original Cameretti correlation should be used.

[pw0908]

Just wanted to give you both an update: I found out what is wrong with all our implementations of ePC-SAFT. Turns out, within the DH theory, the temperature-dependent hard-sphere diameter is used, rather than just sigma. I managed to get this clarified after a few emails with Christoph Held. I am now able to reproduce their results.

[zmeri]

Thank you for sharing this @pw0908. Is this the commit where you fixed it in Clapeyron.jl?
ClapeyronThermo/Clapeyron.jl@bf1c770

I am trying to figure out where the change needs to be made, but I don't see exactly where the temperature dependence comes in. Does the temperature-dependent hard-sphere diameter need to be used for calculating chi and sigma_k in the ion term of ePC-SAFT?

[pw0908]

That is indeed the commit where I fixed it. I allowed the ion diameter to be a kwarg instead where, for ePC-SAFT specifically, we can swap in the temperature-dependent hard-sphere diameter instead.

Indeed, everywhere you would normally see sigma_k in DH, you must now use HSd_k.

[zmeri]

I finally got around to fixing this. Let me know if you see anything I missed when making this fix.

[pw0908]

Did you manage to benchmark to our implementation? I know that we have almost machine precision accuracy with regular PC-SAFT. I’d imagine we’d get quite close with ePC-SAFT as well.

[zmeri]

I tried to compare, but I haven't quite been able to figure out the correct syntax for using the ePCSAFT model in Clapeyron.jl. I think I'm supposed to give the ePCSAFT method an array of solvents and an array of ions, but I keep getting errors. Can you see what I'm doing wrong?

using Clapeyron

model = ePCSAFT(["water"], ["Na+", "Cl-"])
V = 1/55757.07260200306
T = 298.15
p = pressure(model, V, T)
display(p)

[pw0908]

Here is a working example:

using Clapeyron

model = ePCSAFT(["water"], ["sodium", "chloride"])
V = 1/55757.07260200306
T = 298.15
z = [0.90, 0.05, 0.05]
p = pressure(model, V, T, z) # 2.372571785875148e8

[zmeri]

Thanks for the example. For the molar density of toluene I'm getting a very close match. At 320 K and 101325 Pa
Clapeyron.jl: 9033.114359706226 mol/m^3
My PC-SAFT code: 9033.11421467559 mol/m^3

For water, though, I'm getting a bigger difference. Could you help me out with water?

model = PCSAFT(["water"])
z = [1]
T = 274 # K
p = 101325 # Pa
rho_water = molar_density(model, p, T, z)
display(rho_water) # 51961.10187494896 mol/m^3

The IAWPS95 EOS gives 55502.5970532902 mol/m^3 for these conditions, so the value of 51961.10187494896 mol/m^3 from Clapeyron.jl is a -6.4% difference. Am I right that this difference is because a temperature dependent sigma is not used for water? Is there a temperature dependent model for water build-in to Clapeyron.jl?

[pw0908]

Oh, I forgot to mention. Our water model with the temperature-dependent sigma is actually called water08 and we only include the temperature dependence in a variant of PC-SAFT we call pharmaPCSAFT. If you use this, you'll get the correct values:

model = pharmaPCSAFT(["water08"])

molar_density(model,101325,274.) # 55477.149055578426

We did this split just so that we could keep using the default water model in PC-SAFT and allow users to switch when needed (we also noticed that the temperature-dependent variant was only really used for pharmaceutical / biological systems, hence its name).

[zmeri]

Ok yes, now with water08 and pharmaPCSAFT I'm getting comparable results. Here are the checks I've done.

Molar density of water at 274 K and 101325 Pa:
Clapeyron.jl: 55477.149055578426 mol/m^3
My PC-SAFT code: 55476.442745328946 mol/m^3

Molar density of aqueous NaCl at 298.15 K and 101325 Pa:
Clapeyron.jl: 55757.90045018177 mol/m^3
My PC-SAFT code: 55757.07260200306 mol/m^3

Bubble pressure of toluene at 572.6667 K
Clapeyron.jl: 3290657.9543914576 Pa
My PC-SAFT code: 3290656.6556331236 Pa

Bubble pressure of water at 572.6667 K
Clapeyron.jl: 66941.16638384997 Pa
My PC-SAFT code: 66917.70502934992 Pa

And a lot of these results also depend on the algorithm used to solve for things like XA in the association term, density, and bubble pressure.

[pw0908]

Yes these small differences are to be expected. I would recommend trying to compare lower-level functions. For example, we have an a_res function to give the reduced Helmholtz free energy in terms of V, T, z:

model = ePCSAFT(["water"], ["sodium", "chloride"])
V = 1/55757.07260200306
T = 298.15
z = [0.90, 0.05, 0.05]

a_res = Clapeyron.a_res(model, V, T, z) # -9.001262927626707

Comparing these might be better since everything is analytical now. Even the association fraction can be solved exactly.

[zmeri]

Yes, that's a good point. For those same conditions, my PC-SAFT code gives -9.001499164113106 J/mol.

[pw0908]

Ah, I guess this is where we compare parameters 😅. You can get all of ours from the model object:

model = ePCSAFT(["water"], ["sodium", "chloride"])

model.neutralmodel.params.epsilon

model.ionmodel.params.sigma

if I had to guess, are you still using your water model and not the original from Held?

If it’s not, then we might want to check our model / universal constants.

[zmeri]

Ah okay, you want to figure out why the results aren't even closer. I have to admit that since I'm no longer working in the field and it's hard to find time to work on this PC-SAFT code anymore, I was fine with how close the results were before. But, I can check the parameters for you.

Checking the parameters I see that there were some that were slightly different. There were several for which you had more decimal places. I updated the parameters to match when running my code and the result is now closer. My code gave -9.001235762324228 J/mol, when updating the parameters to match yours. Also, for these calculations here I am using the water model for sigma that was given by Held, et al. Here's what I used, in case you want to check.

import numpy as np
from pcsaft import dielc_water, pcsaft_ares

# Aqueous NaCl
# 0: Na+, 1: Cl-, 2: water
x = np.asarray([0.05, 0.05, 0.90])
m = np.asarray([1, 1, 1.2046817736])
s = np.asarray([2.8232, 2.75996, 0.])
e = np.asarray([230.00, 170.00, 353.9449])
volAB = np.asarray([0, 0, 0.04509])
eAB = np.asarray([0, 0, 2425.6714])
k_ij = np.asarray([[0, 0.317, 2.37999],
                    [0.317, 0, -0.25],
                    [2.37999, -0.25, 0]])
z = np.asarray([1., -1., 0.])

t = 298.15 # K
s[2] = 2.7927 + 10.11*np.exp(-0.01775*t) - 1.417*np.exp(-0.01146*t) # temperature dependent segment diameter for water
k_ij[0,2] = -0.007981*t + 2.37999
k_ij[2,0] = -0.007981*t + 2.37999
dielc = dielc_water(t)

den = 55757.07260200306
params = {'m':m, 's':s, 'e':e, 'e_assoc':eAB, 'vol_a':volAB, 'k_ij':k_ij, 'z':z, 'dielc':dielc}
ares_calc = pcsaft_ares(t, den, x, params)
print(ares_calc)

Note that I'm using the version of the code that I have wrapped into a Python package (repo: https://github.com/zmeri/PC-SAFT). It uses the same core C++ functions that are used here in CoolProp, but it's easier to work. For instance, I can more easily change the parameters.

[pw0908]

I didn’t realise you weren’t one of the maintainers of PC-SAFT on CoolProp. I only brought this up cause I know @ibell cares a lot about getting these things as close as possible. This still isn’t quite as close as I would hope, but that might come down to the universal and model constants.

I’ll leave it up to you guys whether or not this is fixed.

[zmeri]

I contributed the PC-SAFT code to CoolProp since I already had the C++ code for PC-SAFT, so that's why I take some of the issues related to PC-SAFT. In recent years though, I've moved on to other projects, so I don't think I contribute enough to really be considered a maintainer.

Thank you for your help with this issue and for taking the time to bring up the error with the ion diameter.

By the way, Clapeyron.jl looks like a pretty slick package. It's impressive to see how many different models and properties are included, and I imagine that attention to getting things as close as possible has paid off. Great work with it!

[pw0908]

No worries! The ion diameter issue was really annoying me for a while (especially considering how not-obvious it is). Glad we could get this sorted!

Thanks a lot for the complement! We've only just started supporting electrolytes; surprisingly, compared to other EoS, the literature regarding these models are quite difficult to reproduce...