CPA example please

[ianhbell]

There is something wrong with my implementation of CPA in teqp (☹️) so I am trying to check my calculations with Clapeyron.jl. The docs could use some units for the model parameters especially a and b: https://clapeyronthermo.github.io/Clapeyron.jl/dev/eos/saft/#Clapeyron.CPA

I am trying to replicate this example with Clapeyron.jl (below code is from with teqp) with the parameters from https://www.sciencedirect.com/science/article/pii/S0378381217300444:

# https://www.sciencedirect.com/science/article/pii/S0378381217300444#tbl8
bi = 182.0037/1e6 # mL/mol -> m^3/mol; 1 mL = 1e-6 m^3
a0i = 10.6863 # Pa m^6/mol^2
c1 = 0.8413
Tc = 643.2 
epsABi = 2679.3736*R # epsABi/R [K], multiply by R to get [J/mol]
beta = 0.0912

pure = {
    "a0i / Pa m^6/mol^2": a0i, "bi / m^3/mol": bi, "c1": c1,
    "Tc / K": Tc, "epsABi / J/mol": epsABi, "betaABi": beta, "class": "2B"
}

j = {"cubic": "SRK", "pures": [pure], "R_gas / J/mol/K": R}
model = teqp.make_model({'kind':'CPA', 'model':j}, validate=False)

rho_molm3 = 5230.088495575221 # m^3/mol
print(get_p_PVT(model=model, T=333, rho=rho_molm3)/1e5, 'bar')

which should be giving me about 1 bar but is giving me 470 bar 🙈

[longemen3000]

@pw0908 should know better than me, but i look up the paper cited in the docs (doi.org/10.1021/ie051305v), and it has the units that our CPA uses:

[ibell]

So far I tried this (to build a CPA model by name and then to copy the parameters into a new one which should avoid unit conversion issues I hope(d) ) to test that I know how to do this, and it didn't go well. I tried so far:

using Clapeyron
m1 = Clapeyron.CPA(["Methanol"])
m2 = Clapeyron.CPA(["Methanol"]; userlocations = (;
                    a = [m1.params.a],
                    Tc = [m1.params.Tc],
                    b = [m1.params.b]),)

which gives the error

KeyError: key "sites" not found

Stacktrace:
 [1] getindex
   @ .\dict.jl:484 [inlined]
 [2] CPA(components::Vector{String}; idealmodel::Type, cubicmodel::Type{RK}, alpha::Type, mixing::Type, activity::Nothing, translation::Type, userlocations::NamedTuple{(:a, :Tc, :b), Tuple{Vector{PairParam{Float64}}, Vector{SingleParam{Float64}}, Vector{PairParam{Float64}}}}, ideal_userlocations::Vector{String}, alpha_userlocations::Vector{String}, activity_userlocations::Vector{String}, mixing_userlocations::Vector{String}, translation_userlocations::Vector{String}, verbose::Bool, assoc_options::AssocOptions)
   @ Clapeyron C:\Users\ian\.julia\packages\Clapeyron\XyeVJ\src\models\SAFT\CPA\CPA.jl:89
 [3] top-level scope
   @ In[2]:1

and I have no idea where/how I should be specifying the sites.

[longemen3000]

Oh?, pass bondvol = nothing and epsilon_assoc = nothing (it is a bug that I'm fixing in the next release)

[pw0908]

Hi Ian,
I briefly gave it a shot in the following example:

julia> model = CPA(["[EMIM][BF4]"]; userlocations=(;
            a = [10.6863/1e1],
            b = [182.0037/1e3],
            c1 = [0.8413],
            Mw = [1.],
            Tc = [643.2],
            Pc = [1e6],
            n_H = [1],
            n_e = [1],
            epsilon_assoc = Dict((("[EMIM][BF4]","H"),("[EMIM][BF4]","e"))=>2679.3736/1e2*8.314),
            bondvol = Dict((("[EMIM][BF4]","H"),("[EMIM][BF4]","e"))=>0.0912*1e3)),
       alpha_userlocations=(;
            c1 = [0.8413],
            Tc = [643.2]))

julia> 1/volume(model,1e5,333.)
100.84175953851668

Our implementation doesn't seem to agree with the author's value (or your unfortunately). I believe I did all the unit conversions correctly (honestly, we should probably just stick to SI when it comes to CPA; I was under the impression users of CPA preferred the units given in the paper).

As for what could be responsible between the various implementations, I have a feels the authors of that paper may have used a different version of CPA (yup, a few variants of it exist in the literature). The values given do seem a bit odd to me as ionic liquids should been highly non-volatile species, with these parameters being a little small to achieve that.

Between our implementations, perhaps we should compare the cubic and association contributions?

[ibell]

Oy - I didn't realize there were multiple versions of CPA. That's not so good. What I implemented matches Python code from my colleagues, but I do wonder about the parameters, it could be well a subtly different model formulation (or bug(s) in teqp is always a possibility). Did you follow exactly the implementation from Kontogeorgis referenced in the docs? When you have a pure fluid, do you specialize to calculate the site fractions without iteration following Huang and Radosz?

Following your example, I tried to copy parameters from one model to another, like so:

m2 = CPA(["Methanol"]; userlocations=(;
    a = m1.params.a.values,
    b = m1.params.b.values,
    c1 = m1.params.c1.values,
    Mw = m1.params.Mw.values,
    Tc = m1.params.Tc.values,
    Pc = [1e6],
    n_H = [1],
    n_e = [1],
    epsilon_assoc = Dict((("Methanol","H"),("Methanol","e"))=>244.76),
    bondvol = Dict((("Methanol","H"),("Methanol","e"))=>15.4)),
    alpha_userlocations=(; c1 = m1.params.c1.values, Tc = m1.params.Tc.values)
    )

but it didn't work, the model can be constructed, but the values obtained from this are very different:

I thought maybe it was the number of sites but that seems identical too.

Note: I couldn't figure out how to index the values for bondvol or epsilon_assoc to copy them over, so I copied them manually.

[longemen3000]

we do perform some transformations on a and b. ideally, we should provide a function like

m2 = CPA(["methanol"],userlocations = to_userlocations(m1))

that allows transforming the cached values into a common form (work is being done, in particular to get/get k and l values, and there is a export_model functionality already available with some models). the inverse transformations of params in the case of CPA are the the following:

m2 = CPA(["Methanol"]; userlocations=(;
    a = m1.params.a.values .* 10, #undo scaling of a
    b = m1.params.b.values .* 1000, #undo scaling of b
    c1 = m1.params.c1.values,
    Mw = m1.params.Mw.values,
    Tc = m1.params.Tc.values,
    Pc = m1.cubicmodel.params.Pc.values,
    n_H = [1],
    n_e = [1],
    epsilon_assoc = Dict((("Methanol","H"),("Methanol","e"))=>244.76),
    bondvol = Dict((("Methanol","H"),("Methanol","e"))=>15.4)),
    )

with that, both models give the same result:

[ianhbell]

I think I agree with @pw0908 - it would be best if CPA coefficients were in SI units. Then it is unambiguous (or at least less ambiguous) what the units are. That is for instance what is done in teqp. I've learned the hard way to only use base SI units everywhere.

Thanks @longemen3000 for the example; I didn't know that the values were rescaled internally.

[longemen3000]

btw, i looked the references and it seems that the paper uses the simplified CPA term (g(η) = 1/(1-1.9η), sCPA in Clapeyron). still cannot reproduce

[ianhbell]

I missed this note in the thread :( See my conclusions at the end.

[ianhbell]

Also, the parameters in Kontogeorgis (what is cited for CPA) are intended for use with sCPA, not CPA (if I am reading properly)

[ianhbell]

@longemen3000 don't sweat it; its now my problem

[pw0908]

Just as a small update / sanity check, I went ahead and re-benchmarked our CPA implementation and can confirm it reproduces the results for pure water as published in the water review paper that Georgios recently published. As such, I think we can safely say that something is going wrong within the implementation of the authors of that paper you shared. As I mentioned earlier, the actual values themselves don't seem like what one might expect for ionic liquids.

[ianhbell]

Can you give some numerical values? I tried to replicate the calculations with teqp and I get something quite different.

In Julia I did:

using Clapeyron
using Unitful
water = Clapeyron.CPA(["Water"])
volume(water, 1e5u"Pa", 303.15u"K")

yielding

1.7915123921401366e-5 m³

Then I tried to do the same calculation in teqp, taking the Clapeyron value of density, with the parameters taken from
https://doi.org/10.1016/S0378-3812(99)00060-6 (the ones in Clapeyron seem to be truncated to match the table in Kontogeorgis: https://pubs.acs.org/doi/epdf/10.1021/ie051305v)

import teqp, numpy as np
R = 8.31446261815324
water = {
    "a0i / Pa m^6/mol^2": 0.12277 , "bi / m^3/mol": 0.000014515, "c1": 0.67359,
    "Tc / K": 647.096, "epsABi / J/mol": 16655.0, "betaABi": 0.0692, "class": "4C"
}
j = {"cubic": "SRK", "pures": [water], "R_gas / J/mol/K": R}

model = teqp.make_model({'kind':'CPA', 'model':j}, validate=False)
rho = 1/1.7915123921401366e-5
T = 303.15
p = rho*R*T*(1+model.get_Ar01(T, rho, np.array([1.0])))
p

yielding

10923.734492239215

which is the back-calculated pressure in Pa, which is quite different. I guess it is time to debug into the codes to see where the divergence happens.

Is there a way to obtain the association fractions somehow from Clapeyron? That would be the first obvious place to look for bugs in my code.

[ianhbell]

P.S. The rounding on bi seems to make an enormous difference

[ianhbell]

Ok, I understand what is going on here I think. Everything checks out, but there is a key thing that needs to be revised in Clapeyron: the radial dist function needs to be allowed to be user-specified. I tried the two g functions, and if I switch back to the original version of g (like what is implemented in Clapeyron), I get the pressure to be 100 kPa. I had implemented the g function from sCPA as given in the Kontogeorgis paper. Changing only that makes an enormous difference: a factor of 70 for this liquid state. These are the two functions:

The second is what was implemented in teqp, the first what is implemented in Clapeyron.

[ianhbell]

If I understand their paper, they seem to suggest they used the Eq. 5a formulation

[longemen3000]

I will have to check the parameters, because, while the CPA parameters are similar, they are not the same as the ones appearing in the Kontogeorgis paper. we already have sCPA as an EoS, but i saw in the teqp code that there are 2 additional rdf functions (that would be useful to have here too)

[ianhbell]

I want to highlight here the importance of the right radial distribution function and parameters. I tried two values for b that differ only a little bit (3 significant digits versus 5), and the two radial distribution functions:

import teqp, numpy as np
R = 8.31446261815324

for bi in [0.000014515, 0.0000145]:
    for radial_dist in ['KG', 'CS']:
        water = {
            "a0i / Pa m^6/mol^2": 0.12277 , "bi / m^3/mol": bi, "c1": 0.6736,
            "Tc / K": 647.13, "epsABi / J/mol": 16655.0, "betaABi": 0.0692, "class": "4C"
        }
        j = {"cubic": "SRK", "pures": [water], "R_gas / J/mol/K": R, "radial_dist": radial_dist}

        model = teqp.make_model({'kind':'CPA', 'model':j}, validate=False)
        rho = 1/1.7915123921401366e-5
        T = 303.15
        p = rho*R*T*(1+model.get_Ar01(T, rho, np.array([1.0])))
        print(f'{radial_dist}: {p/1e5} bar @ bi= {bi} m^3/mol')

yielding

KG: 109.15910698115533 bar @ bi= 1.4515e-05 m^3/mol
CS: 32.47504686042741 bar @ bi= 1.4515e-05 m^3/mol
KG: 78.17651421286142 bar @ bi= 1.45e-05 m^3/mol
CS: 1.0000000000114575 bar @ bi= 1.45e-05 m^3/mol

So the differences are huge!

[longemen3000]

i'm working on standardizing CPA, and implementing a radial_dist::Symbol argument to select the RDF. i saw in your code that you implement the OT argument, isn't that the same that KG ? (0.475*4 = 1.9). it is relevant in the context of Cubic + Chain + Assoc (CPCA)

[ianhbell]

So it is... So the OT option is duplicative, I'll drop it.

[longemen3000]

on the CPA_SI branch: