Induced-Association Interactions in SAFT

[1578jason]

Hi, I'm trying to use the SAFT EOS to do some calculations on "induced solvation", where only one compound is self-associating but can interact with another compound. I tried to do the calculations as described in the following two articles, but was confused about how to define these systems in Clapeyron, so I come to you for help:
Article 1 https://doi.org/10.1021/jp072640v
Article 2 https://doi.org/10.1021/ie050976q
In the first article, the authors make the following assumption:
• The association energy parameter of the non-self-associating component is set to zero.
• The association volume parameter of this component is assumed to be equal to the value of the associating component in the mixture.
Is it represented in Clapeyron as something like this? (Take ethanol and chloroform as example)

model = PCSAFT(["ethanol","chloroform"];userlocations=(;
Mw = [46.069,119.37],
segment = [2.3827,2.5038],
sigma = [3.1771,3.4709],
epsilon = [198.24,271.625],
n_H=[1,1],
n_e=[1,1],
epsilon_assoc = Dict(
    (("ethanol","e"),("ethanol","H")) =>2653.4,
    (("ethanol","H"),("chloroform","e")) =>1326.7,
    (("ethanol","e"),("chloroform","H")) =>1326.7,
    (("chloroform","H"),("chloroform","e")) =>0.
    ),
bondvol = Dict(
    (("ethanol","e"),("ethanol","H")) => 0.032384,
    (("ethanol","H"),("chloroform","e")) => 0.032384,
    (("ethanol","e"),("chloroform","H")) => 0.032384,
    (("chloroform","H"),("chloroform","e")) => 0.032384
)
))

In the second paper, the authors set the cross-association energy to half of the self-associating component, and take the cross-association volume and kij as two adjustable parameters, whether it is expressed in the following form?(Specify only one pair of cross-association parameters and regress cross-association volume )

model = PCSAFT(["ethanol","chloroform"];userlocations=(;
Mw = [46.069,119.37],
segment = [2.3827,2.5038],
sigma = [3.1771,3.4709],
epsilon = [198.24,271.625],
n_H=[1,1],
n_e=[1,0],
epsilon_assoc = Dict(
    (("ethanol","e"),("ethanol","H")) =>2653.4,
    (("ethanol","e"),("chloroform","H")) =>1326.7,
    ),
bondvol = Dict(
    (("ethanol","e"),("ethanol","H")) => 0.032384,
    (("ethanol","e"),("chloroform","H")) => 0.032384,
)
))

[pw0908]

Hi! A few things to note:

1. The first paper you refer to uses PCP-SAFT which is available in Clapeyron (you'll need to specify the dipole moment in those cases).
2. It appears that the approach in both papers, as far as the induced association is concerned, is identical, where both use the Wong-Sandler combining rule for cross-association.
3. In the first paper, in the case of Chloroform, the visual they've included is a bit misleading. Chloroform is not a proton donor (that C-H bond is not polar enough to induce hydrogen bonding) but it is an electron donor (the C-Cl bond is very polar and can induce hydrogen bonding with the O-H on the alcohols). As there are three of these, you can account for them as e sites.

With those two points in mind, the 'correct' way to set-up the model would be as follows:

julia> model = PCPSAFT(["ethanol","chloroform"];userlocations=(;
       Mw = [46.069,119.37],
       segment = [2.3827,2.5038],
       sigma = [3.1771,3.4709],
       epsilon = [198.24,271.625],
       dipole = [0.,1.04],
       n_H=[1,3],
       n_e=[1,0],
       epsilon_assoc = Dict(
           (("ethanol","e"),("ethanol","H")) =>2653.4,
           (("ethanol","e"),("chloroform","H")) =>1326.7,
           ),
       bondvol = Dict(
           (("ethanol","e"),("ethanol","H")) => 0.032384,
           (("ethanol","e"),("chloroform","H")) => 0.032384,
       )
       ))
PPCSAFT{BasicIdeal} with 2 components:
 "ethanol"
 "chloroform"
Contains parameters: Mw, segment, sigma, epsilon, dipole, dipole2, epsilon_assoc, bondvol

julia> bubble_temperature(model,1e5,[0.8,0.2])
(342.47489071408995, 6.664468359724602e-5, 0.027564634355155437, [0.5689549253463148, 0.43104507465368525])

The values above appear to correctly reproduce the results from the paper although I haven't rigorously tested this.

I hope this helps!

[1578jason]

Hi Pierre,
Thanks for your quick reply.
In your suggestion, I note that chloroform is treated as an electron donor substance, Then it seems that it should be changed to?

n_H=[1,0],
n_e=[1,3],

so in this case ,the cross-association exists only between the proton sites of ethanol and the electron sites of chloroform. However, I found another paper on the early (10.1016 / j.fluid. 2007.04.015).In this paper, 1,1,1,2,3,3,3-heptafluoropropane is considered to have two association sites (one proton site and one electron site) just like ethanol. Then, using the assumption in the first paper, I reconstructed the model in Clapeyron to accurately reproduce the results in the article.

model = PCSAFT(["ethanol","HFC-227ea"];userlocations=(;
Mw = [46.069,119.37],
segment = [2.3827,3.5190 ],
sigma = [3.1771,3.3165 ],
epsilon = [198.24 182.573;
        182.573 164.18],
n_H=[1,1],
n_e=[1,1],
epsilon_assoc = Dict(
    (("ethanol","e"),("ethanol","H")) =>2653.4,
    (("ethanol","e"),("HFC-227ea","H")) =>1326.7,
    (("ethanol","H"),("HFC-227ea","e")) =>1326.7,
    (("HFC-227ea","e"),("HFC-227ea","H")) =>0.,
    ),
bondvol = Dict(
    (("ethanol","e"),("ethanol","H")) => 0.032384,
    (("ethanol","e"),("HFC-227ea","H")) => 0.032384,
    (("ethanol","H"),("HFC-227ea","e")) => 0.032384,
    (("HFC-227ea","e"),("HFC-227ea","H")) => 0.032384,
)
))
x=[0.,0.0604,0.1228,0.3314,0.5219,0.7260,0.8547,0.9440,1]
for i=1:9
bub=bubble_pressure(model,343.13,[x[i],1-x[i]])
println(bub[1])
end
1.485265522487281e6
1.4324384088390123e6
1.3864812966939074e6
1.2691594695010856e6
1.139143438538306e6
874437.2146160462
581260.7514751055
294315.16363260965
71805.77594729989

It seems that this definition model is also feasible?

Best regards,
jason

[pw0908]

Hi Jason,

Good catch on my typo! It shouldn't have made a difference in the calculations.

Looking at this new paper you've sited, it makes mention of "HFC-227ea does not self-associate but act as proton donor in the mixture with ethanol", nothing about it acting as an electron donor. At the end of the day, since this association parameters are identical, it shouldn't have a large impact on the actual phase equilibrium. The more-correct form for me would have been 2 proton sites and 0 electron sites, but if this reproduces the results from the paper, then I wouldn't worry about this too much.

You seem to have gotten a hang of it!

Best,

Pierre

[1578jason]

Hi Pierre,
I'm using the parameter estimation tool to fit several sets of kij, but I encountered the following error when modeling the methane-carbon dioxide system. Can you help me take a look?

using Clapeyron, BlackBoxOptim,Metaheuristics

model = SAFTVRMie(["methane","carbon dioxide"]; userlocations=(;
Mw = [16.04,44.01], 
segment = [1,1.5],
sigma = [3.737,3.1916], 
epsilon = [152.58,231.88], 
lambda_a = [6.,5.1646],
lambda_r = [12.504,27.557],
epsilon_assoc =nothing,
bondvol = nothing
))

toestimate = [
    Dict(
        :param => :epsilon,
        :indices => (1,2),
        :symmetric => true,
        :lower => 150.,
        :upper => 240.,
        :guess => 186.
    )
];

function bubble_point(model::EoSModel,T,x)
    bub = bubble_temperature(model,T,[x,1-x])
    return bub[1], bub[4][1]
end

estimator,objective,initial,upper,lower = Estimation(model,toestimate,["data/methane+carbon dioxide.csv"]);

function logger(st)
    if st.iteration % 10 == 0
        sol = st.best_sol    
        it = st.iteration
        println("Iteration: ",st.iteration)
        println("Best solution: ",st.best_sol)
    end
end

bounds = [lower upper]'
result = optimize(objective, bounds, ECA(),logger=logger)
params = minimizer(result)

Errors are as follows:

ERROR: MethodError: no method matching /(::Vector{Missing}, ::Missing)
Closest candidates are:
  /(::AbstractArray, ::Unitful.Units) at C:\Users\.julia\packages\Unitful\J4AJj\src\quantities.jl:46
  /(::AbstractVecOrMat, ::PDMats.ScalMat) at C:\Users\.julia\packages\PDMats\bzppG\src\scalmat.jl:55
  /(::AbstractVecOrMat, ::PDMats.PDMat) at C:\Users\.julia\packages\PDMats\bzppG\src\pdmat.jl:65

data.zip

Best regards,

Jason

[pw0908]

Hi Jason,

Having taken a look at your code, I've got a few comments:

1. The issue you were having actually has now to do with the code itself. If you open your CSV file in VSCode, you'll see that there are about 100 blank rows. These are included as part of your data set and, since they are empty, lead to the error you were seeing here. The simple fix is to delete those rows.
2. I took a look at how you were fitting your data. You were using pxy data (so constant temperature) but fit to the temperature and vapour composition. I'd recommend fixing the temperature instead and fitting to the pressure and vapour fraction. This is probably more-representative of the experimental data, and will be a lot faster. I've made these changes to your script attached bellow.
3. Another issue I encountered when trying to run this fitting is that your data includes points very close to the critical point of the mixture. It is generally recommended, when fitting SAFT-type equations, to fit away from the critical point. This is mainly because, for most SAFT equations, the underlying theory inherently fails / worsens near the critical point. Removing some of this data might improve the stability of the calculations as, right now, these are very slow. I had to switch the method being used (FugBubblePressure()) to improve the stability.
4. Finally, and perhaps a little bit anti-climatically, you'll be happy to know that the carbon dioxide-methane unlike parameters have already been fitted to the very same experimental data you are using in the following paper: https://doi.org/10.1016/j.fluid.2015.12.041 . The pure-fluid parameters they are using are not the same as yours, but, if you're just interested in modelling this phase equilibrium specifically, I'd recommend just using these parameters.

I've included the scripts I've used below.
test.zip

Hope it helps!

• Pierre

[1578jason]

Hi! I recently had an idea to compare the computational efficiency of models in Clapeyron, such as sPC-SAFT versus the original PC-SAFT, or even simple cubic equations versus complex SAFT models (low computational efficiency is one of the reasons why SAFT is difficult to generalize). I plan to start with the basic binary phase equilibrium comparision and focus on evaluating the running time required for the models to track the phase envelope.However, I am unsure if the code in the notebook examples provides a fair comparison? Are the conditions, other than the model itself, such as the algorithm, consistent across the comparisons? Or any other better implementations in Clapeyron that could be considered? Any suggestions and guidance would be greatly appreciated!

[pw0908]

Hi Jason! If you're going to do any benchmarking, I recommend using BenchmarkTools.jl which will efficiently, and fairly, sample the computational times. Here is an example:

using Clapeyron, BenchmarkTools
model = PR(["water","ethanol"])

@benchmark bubble_pressure(model,298.15,[0.5,0.5])
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  60.959 μs …   3.645 ms  ┊ GC (min … max): 0.00% … 96.33%
 Time  (median):     64.625 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   72.792 μs ± 158.701 μs  ┊ GC (mean ± σ):  9.78% ±  4.40%

            ▁ ▃▅█▄▄▂                                            
  ▁▂▃▃▄▄▄▄▆▇████████▇▇▅▅▃▄▄▄▃▄▃▃▃▃▃▂▃▃▂▂▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▃
  61 μs           Histogram: frequency by time         75.3 μs <

 Memory estimate: 95.34 KiB, allocs estimate: 917.

As far as binary phase equilibrium calculations go, whether its cubics or SAFT, they go through the same algorithms (@longemen3000 can correct me on this). With that in mind though, if you're trying to highlight the advantages of cubics compared to SAFT equations, wouldn't it be better to look at properties like bulk volume calculations? Cubics have an analytical solution whereas SAFT equations need to solved for iteratively:

# PR
model = PR(["water","ethanol"])
@benchmark volume(model,1e5,298.15,[0.5,0.5])
BenchmarkTools.Trial: 10000 samples with 72 evaluations.
 Range (min … max):  840.847 ns … 102.039 μs  ┊ GC (min … max): 0.00% … 98.48%
 Time  (median):     886.583 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   938.731 ns ±   1.907 μs  ┊ GC (mean ± σ):  4.66% ±  2.40%

           ▁▃▅▄█▄▃▅▂▁▃  ▁                                        
  ▁▁▁▂▂▄▄▅███████████████▇█▇▆█▆▆▆▄▅▅▄▃▃▃▃▂▃▂▂▂▂▂▂▂▁▂▂▁▂▁▂▁▁▁▁▁▁ ▃
  841 ns           Histogram: frequency by time          995 ns <

 Memory estimate: 400 bytes, allocs estimate: 5.

# PC-SAFT
model = PCSAFT(["water","ethanol"])
@benchmark volume(model,1e5,298.15,[0.5,0.5])
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  126.250 μs …  2.102 ms  ┊ GC (min … max): 0.00% … 90.52%
 Time  (median):     142.792 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   147.146 μs ± 48.532 μs  ┊ GC (mean ± σ):  0.63% ±  2.02%

           ▄▇█▆█▇▄▂▁▁▂▂▃▁                                       
  ▂▁▁▁▂▂▃▄▇███████████████▆▅▄▄▄▃▃▃▃▃▂▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂ ▄
  126 μs          Histogram: frequency by time          191 μs <

 Memory estimate: 25.39 KiB, allocs estimate: 116.

Thats a 3-order of magnitude difference!

[pw0908]

And actually, you could also consider pure-component saturation conditions where cubics and PC-SAFT have Super-Ancillary equations. These provide speed-ups in estimating saturation points (for PC-SAFT):

model = PCSAFT(["hexane"])

# Iterative method
@benchmark saturation_pressure(model,298.15)
BenchmarkTools.Trial: 10000 samples with 8 evaluations.
 Range (min … max):  3.760 μs … 138.859 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     3.859 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   3.962 μs ±   1.582 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

   ▃▆███▇▇▇▆▅▃▁   ▁      ▂▃▃▃▃▃▃▂▁▁                           ▂
  ███████████████████▇▆▇████████████▇▇▇▇▆▆▆▆▆▆▅▂▄▄▅▅▅▃▅▅▅▅▅▅▅ █
  3.76 μs      Histogram: log(frequency) by time      4.69 μs <

 Memory estimate: 32 bytes, allocs estimate: 1.

# Using superancillary
@benchmark saturation_pressure(model,298.15,SuperAncSaturation())
BenchmarkTools.Trial: 10000 samples with 10 evaluations.
 Range (min … max):  1.663 μs …  25.087 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     1.717 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   1.725 μs ± 270.360 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

           ▁▃ ▅▇▇█▇ ▅▃▁                                        
  ▂▁▂▂▂▁▃▄▆██▁█████▁███▇▅▁▄▃▂▂▂▁▂▂▂▂▂▁▂▂▂▁▂▁▂▁▂▂▂▁▂▂▂▂▂▁▂▂▂▂▂ ▃
  1.66 μs         Histogram: frequency by time        1.87 μs <

 Memory estimate: 32 bytes, allocs estimate: 1.

[1578jason]

Hi Pierre! Thanks for your quick reply! I actually want to look at the correlation between the accuracy of the model and the speed of the calculation (the increased accuracy of the SAFT and the increased cost of the calculation). As has been done in much of the literature, accuracy is often quantified by calculating the deviation of one or more sets of model-calculated values from experimental values. Therefore, I might prefer to get a tracked set of phase boundaries to measure the computational cost of the model. So is it possible to get a fair comparison with the examples in the notebook?

model1=PR(["ethanol","water"])
model2=PCSAFT(["ethanol","water"])
T =303.15
x = range(1e-5,1-1e-5,length=100)
X = Clapeyron.FractionVector.(x)
@benchmark begin
	bub = [] 
        v0=[]
	for j ∈ 1:100
        if j==1
            append!(bub, [bubble_pressure(model2,T,X[j])])
            v0 = [log10(bub[j][2]),log10(bub[j][3]),bub[j][4][1],bub[j][4][2]]
        else
            append!(bub, [bubble_pressure(model2,T,X[j];v0=v0)])
            v0 = [log10(bub[j][2]),log10(bub[j][3]),bub[j][4][1],bub[j][4][2]]
        end
    end
end

[pw0908]

I would advise against benchmarking for-looped operations. The computational cost of an eos is usually defined as the time take for a single operation not to trace a full phase diagram. This is mainly because it is arbitrary how many points one can use to trace the diagram (the above example uses 100, but it could have easily been 50 or 1000). As such, all you really want to do it benchmark the computational cost of a single bubble_pressure calculation.

[longemen3000]

I agree with Pierre in that regard, if you want to compare efficiency, you need to evaluate bulk properties. equilibrium properties can be measured in terms on how much eos evaluations you require to reach an equilibrium state. For example a bubble pressure algorithm my require less iterations, but more EoS evaluations because of repeated volume calculations (FugBubblePressure vs ChemPitBubblePressure), but if you have an specialized volume implementation (like cubics), then the cost model flips, and you want to perform as much volume evaluations as possible.

With regards to supperancillaries, it seems like cheating, but in my opinion, if they are available, they should be used as much as possible. The only reason we don't use those by default is because of size (the EoSSuperancillaries.jl PCSAFT superanc parameters weight around 3mb, and the REFPROP superanc parameters weight around 12 mb).

Also, due to how phase envelopes are mathematically, you will get convergence problems (and more EoS evaluations) near the critical mixture point with any EoS, that is just how critical points work.

Finally answering @pw0908 question, there are some differences on binary equilibrium, but just structural changes (activity models need to be solved with the activity solver). the main changes between cubics and SAFT are:

• cubics have a direct implementation of crit_pure (for obvious reasons), SAFTs require iterating (we try to minimize critical calculations as much as possible)
• cubics have a volume implementation that only requires 1 eos evaluation, SAFTs require iterations.
• due to the first two reasons, the initial points for saturation pressure for cubics are faster (around 1 eos evaluation for cubics vs at least 1 volume iteration + some additional eos evaluations for other EoS). This last point does not matter if you are using superancillaries.

[1578jason]

Ok, I have no more questions, thanks again for your advice and help!

[1578jason]

Hi! I want to repeat the work of this paper(10.1016/j.fluid.2020.112646), so I wonder if there is a quick way to replace the universal constants to creat a new model. Thanks!

[pw0908]

This is actually pretty straight-forward. You'll just need to define a new EoS (LKPCSAFT) with the appropriate functions modified:

using Clapeyron
import Clapeyron: SA, PCSAFTModel, IdealModel, PCSAFTParam, @newmodel

abstract type LKPCSAFTModel <: PCSAFTModel end
@newmodel LKPCSAFT LKPCSAFTModel PCSAFTParam{T}

export LKPCSAFT

function I(model::LKPCSAFTModel, V, T, z, n, _data=@f(data))
    _,_,_,_,η,m̄ = _data
    if n == 1
        corr = LKPCSAFTconsts.corr1
    elseif n == 2
        corr = LKPCSAFTconsts.corr2
    end
    res = zero(η)
    m̄1 = (m̄-1)/m̄
    m̄2 = (m̄-1)/m̄*(m̄-2)/m̄
    @inbounds for i ∈ 1:length(corr)
        ii = i-1 
        corr1,corr2,corr3 = corr[i]
        ki = corr1 + m̄1*corr2 + m̄2*corr3
        res +=ki*η^ii
    end
    return res
end

const LKPCSAFTconsts = (
    corr1 =
    SA[(0.9105631445,-0.3084016918, -0.0906148351),
    (0.6361281449, 0.1860531159, 0.4527842806),
    (2.6861347891, -2.5030047259, 0.5962700728),
    (-26.547362491, 21.419793629, -1.7241829131),
    (97.759208784, -65.255885330, -4.1302112531),
    (-159.59154087, 83.318680481, 13.776631870),
    (91.297774084, -33.746922930, -8.6728470368)], # Replace these constants

    corr2 =
    SA[(0.7240946941, -0.5755498075, 0.0976883116),
    (2.2382791861, 0.6995095521, -0.2557574982),
    (-4.0025849485, 3.8925673390, -9.1558561530),
    (-21.003576815, -17.215471648, 20.642075974),
    (26.855641363, 192.67226447, -38.804430052),
    (206.55133841, -161.82646165, 93.626774077),
    (-355.60235612, -165.20769346, -29.666905585)] # Replace these constants
)

You'll need to define parameters manually though to use it:

model = LKPCSAFT(["a1"],userlocations = (;
               Mw = [15.],
               epsilon = [200.;;],
               sigma = [3.;;],
               segment = [1.],
               n_H = [0],
               n_e = [0],
               epsilon_assoc = nothing,
               bondvol = nothing))

You can also do something similar with the density-dependent hard-sphere diameter.

[1578jason]

Thanks for your quick reply! I compare the newly defined model with the original model and find that the calculation results using the same parameters are the same. Could you help me have a look?
case.zip

[longemen3000]

oh, that is just a typo in the definition of I. supposing you are working on the REPL, then it should be:

function Clapeyron.I(model::LKPCSAFTModel, V, T, z, n, _data=@f(Clapeyron.data))

with that, all three models give different results:

julia> bubble_pressure(model1,510.0,[0.5,0.5])
(7.6618698226611875e6, 0.000206589101703951, 0.00040928651781706123, [0.8374865558859176, 0.16251344411408242])

julia> bubble_pressure(model2,510.0,[0.5,0.5])
(8.191676980779865e6, 0.00020567417060249771, 0.0003846972543888197, [0.8385613999174941, 0.16143860008250588])

julia> bubble_pressure(model3,510.0,[0.5,0.5])
(7.935139589073927e6, 0.0002041729768480577, 0.0004016998955630544, [0.8453174645442335, 0.1546825354557665])