zmeri / PC-SAFT

Functions implementing the PC-SAFT equation of state, including association, electrolyte and dipole terms
GNU General Public License v3.0
46 stars 18 forks source link

Ionic liquids density from ePC-SAFT not matching literature #107

Closed DiegoTMelfi closed 1 year ago

DiegoTMelfi commented 1 year ago

Dear,

I was trying to reproduce some results from Ji, Held, and Sadowski (Fluid Phase Equilibria 335 (2012) 64–73). For the density of ionic liquids. From the reference, the reported IL parameters should closely match the experimental data. However, the ePC-SAFT implementation is consistently giving me densities higher than the experimental values. For instance for [C2min][BF4], I have this piece of code:

# [C2min][BF4] - Parameters from Ji, Held, and Sadowski (Fluid Phase Equilibria 335 (2012) 64–73)
# 0 = [C2min](+), 1 = [BF4](-)
MW = np.asarray([111.17, 86.81]) #g/mol
x = np.asarray([0.5, 0.5])
m = np.asarray([1.4872, 3.8227])
s = np.asarray([3.5926, 3.5088])
e = np.asarray([206.49, 496.12])
volAB = np.asarray([0, 0])
eAB = np.asarray([0, 0])
k_ij = np.asarray([[0, 0],
                   [0, 0]])
z = np.asarray([1, -1])
dielc = 1

pyargs = {'m':m, 's':s, 'e':e, 'e_assoc':eAB, 'vol_a':volAB, 'k_ij':k_ij, 'z':z, 'dielc':dielc}

T = 313.2 # K
P = 1e5 # Pa

mol_density = pcsaft_den(T, P, x, pyargs, phase='liq') #mol/m3
aveMW = sum(x*MW) #g/mol
mass_density = mol_density*aveMW*1e-3 #kg/mol

ref = 1269 #kg/m3 from Journal of Chemical & Engineering Data, Vol. 54, No. 1, 2009

print('\n##########  Test for [C2min][BF4]  ##########')
print('----- Density at',T,'K and ',P,' Pa')
print('    Reference:', ref, 'mol m^-3')
print('    PC-SAFT:', mass_density, 'kg/m^-3')
print('    Relative deviation:', (mass_density-ref)/ref*100, '%')

That gives me the following output:

########## Test for [C2min][BF4] ########## ----- Density at 313.2 K and 100000.0 Pa Reference: 1269 mol m^-3 PC-SAFT: 1423.9168766682528 kg/m^-3 Relative deviation: 12.207791699625913 %

It seems like the ion contribution is a little too high. But maybe I am just missing something.

Thank you in advance for your attention! Diego

defencedog commented 1 year ago

@DiegoTMelfi Can u reproduce & check your densities using feos

DiegoTMelfi commented 1 year ago

@DiegoTMelfi Can u reproduce & check your densities using feos

@defencedog As of now, I believe feos doesn't have the ion term implemented. So I cannot use it to reproduce and check the results.

zmeri commented 1 year ago

Hi @DiegoTMelfi, I don't have the time right now to look into this myself and I'm not sure exactly what the problem could be. One potential cause could be the density solver. PC-SAFT often has more than 3 density roots and it can be difficult to find the correct one (see Privat et al. 2010). I have seen before that sometimes the solver will find the wrong root, especially with these kinds of ionic mixtures. To see if that is the case, you could create a small Python function to calculate the pressure for a range of rho values. Then you can calculate the difference from the actual pressure value in your example and see where the roots appear.

Unfortunately, I haven't found a good way to check the ion contribution values because there isn't another available source code for the ion or polar terms with which to compare. You could talk with Christoph Held and ask him if he could share the internal values for that specific point you're testing. Maybe he would also be willing to compare his code to the code here in this repo. I've talked with Christoph before and he has not been willing to share his code openly, which is unfortunate, especially since sharing research is part of the ethos of academia. However, he was at least willing to give results from his code for a single point so that I could compare them to my own code here. For that one point, my results matched his, if I remember correctly. There may still be some problem though that I never noticed.

DiegoTMelfi commented 1 year ago

Hi All,

I contacted Dr. Held and he was very helpful and helped me find the problem.

For ionic liquids in Ji, Held, and Sadowski (Fluid Phase Equilibria 335 (2012) 64–73), the dispersion contribution is calculated as if the ionic liquids were any other molecule. So the d(T) and e_ij are computed as if they were not ions. Since in the near future I won't be working with small ions, I just commented the lines with conditionals that would change d and e_ij for ions.

for example:

//    if (!cppargs.z.empty()) {
//     for (int i = 0; i < ncomp; i++) {
//         if (cppargs.z[i] != 0) {
//             d[i] = cppargs.s[i]*(1-0.12); // for ions the diameter is assumed to be temperature independent (see Held et al. 2014)
//        }
//    }
// }

and:

//            if (!cppargs.z.empty()) {
//                if (cppargs.z[i]*cppargs.z[j] <= 0) { // for two cations or two anions e_ij is kept at zero to avoid dispersion between // like ions (see Held et al. 2014)

Those lines show up a couple of times, so it is important to comment all of them.

It is not the most elegant solution and it restricts the applicability of the code, so I am not putting in a pull request for this. But if someone faces the same issue in the future, this is an easy fix.

Thank you! Diego

zmeri commented 1 year ago

Thank you for sharing this!

On Wed, Nov 16, 2022, 00:16 DiegoTMelfi @.***> wrote:

Hi All,

I contacted Dr. Held and he was very helpful and helped me find the problem.

For ionic liquids in Ji, Held, and Sadowski (Fluid Phase Equilibria 335 (2012) 64–73), the dispersion contribution is calculated as if the ionic liquids were any other molecule. So the d(T) and e_ij are computed as if they were not ions. Since in the near future I won't be working with small ions, I just commented the lines with conditionals that would change d and e_ij for ions.

for example:

// if (!cppargs.z.empty()) {

// for (int i = 0; i < ncomp; i++) {

// if (cppargs.z[i] != 0) {

// d[i] = cppargs.s[i]*(1-0.12); // for ions the diameter is assumed to be temperature independent (see Held et al. 2014)

// }

// }

// }

and:

// if (!cppargs.z.empty()) {

// if (cppargs.z[i]*cppargs.z[j] <= 0) { // for two cations or two anions e_ij is kept at zero to avoid dispersion between // like ions (see Held et al. 2014)

Those lines show up a couple of times, so it is important to comment all of them.

It is not the most elegant solution and it restricts the applicability of the code, so I am not putting in a pull request for this. But if someone faces the same issue in the future, this is an easy fix.

Thank you! Diego

— Reply to this email directly, view it on GitHub https://github.com/zmeri/PC-SAFT/issues/107#issuecomment-1315932038, or unsubscribe https://github.com/notifications/unsubscribe-auth/AFKXTOFMDR5U44RJ5I5BFCLWIQDSFANCNFSM6AAAAAAR5FUBN4 . You are receiving this because you commented.Message ID: @.***>