bubble_pressure with SAFTgammaMie

[denbigh]

First things first - thank you so much for your wonderful work on this package! I am currently doing research on transport phenomena in high pressure ternary mixtures of CO2/water/cosolvent[methanol/acetonitrile/etc.] and the existence of a simple to use group contribution SAFT implementation in Julia, for prediction of fugacity coefficients, and hopefully for analysis of LLE, is a godsend. I have some interest/background/time to work on multicomponent flash calculations; if there are any algorithms which you are looking to implement, I'd be open to lending a hand.

I also have a brief technical question. I have been playing around with Clapeyron, and I have noticed a potential issue with bubble_pressure calculations when using SAFTgammaMie.

The following code for calculating the bubble pressure works well using PCSAFT:

system = PCSAFT(["methanol","water"])
p = 1e5
T = 383.15
T2 = 443.15
z = [0.5,0.5]
z_LLE = [0.27,0.73]

Clapeyron.bubble_pressure(system,T,z)

It also works well with SAFTVRMie. However, with SAFTgammaMie:

system = SAFTgammaMie(["methanol","water"])
p = 1e5
T = 383.15
T2 = 443.15
z = [0.5,0.5]
z_LLE = [0.27,0.73]

Clapeyron.bubble_pressure(system,T,z)

I get an error:

DimensionMismatch("first array has length 2 which does not match the length of the second, 1.")

I believe this may be traced back to lines 12 and 13 of bubble_pressure.jl, which require splitting the mixture model into pure models:

pure = split_model(model)
crit = crit_pure.(pure)

I noticed that split_model appears to split SAFTgammaMie mixtures into some kind of composite object with type:

SAFTgammaMie{BasicIdeal, SAFTVRMie{BasicIdeal}}

whereas if you use split_model on, say, SAFTVRMie, you get objects of a simpler type:

SAFTVRMie{BasicIdeal}

Do you know if this is a bug in either split_model or in the SAFTgammaMie implementation? It would be useful to do bubble_pressure calculations using SAFTgammaMie.

[longemen3000]

Hi!
i confirm the bug in SAFTgammaMie. The error is not in split_model but in crit_pure, and that is linked at the association energy evaluation (site indices mismatch). i'm gonna try to fix it as soon as possible.
the main culprit it seems that the implementation of the "site squasher" used by SAFTgammaMie generates some indices out of bounds
i'm also looking for good algorithms for multicomponent flash, if you want to contribute with any multiphase algorithm,feel free to ask any questions about the current implementations or make contributions.

[longemen3000]

i think i fixed it. i'm waiting for the tests to pass and then release a new version as soon as possible

[denbigh]

Thanks for this! I have run into similar issues with, e.g., gibbs_free_energy, but I will wait to see if the update fixes these too.

I just want to re-iterate how wonderful Clapeyron.jl is! Even just for using Cubic EOS's, it's a wonderfully simple & efficient interface. I plan to use it for a long time to come.

Regarding multicomponent flash algorithms: to understand the interface, and to answer some of my own research questions, I have now implemented a simple multicomponent, multiphase, isothermal-isobaric flash algorithm, which uses differential evolution (DE) algorithms from BlackBoxOptim.jl to find the global minima of the Gibbs Free Energy, calculated using Clapeyron.jl (c.f. 10.1142/9789813207523_0007). I have tested it on several 2- and 3-species cases, and am yet to see it fail. DE flash calculations will never be as fast as specialized codes, but I have found them to be extremely robust, and also very simple conceptually. The algorithm takes around 10000 function evaluations to converge, which takes about 0.1-5 seconds depending on the underlying EOS. The one downside is the user must specify the number of phases - this is not worked out automatically at present. If you specify too many phases, the code predicts multiple phases with equal composition, and if you specify too few, the predictions are unphysical. You can iteratively solve the code, increasing the number of phases by 1 each time, until the Gibbs Free Energy stops decreasing; this process could probably be automated if desired.

Would there be interest in incorporating a DE flash algorithm into Clapeyron.jl? I'd be happy to formulate it as a pull request if so.

[pw0908]

Hi there!

Great to hear that you plan to use Clapeyron in your future work! I hope these fixes will help you with these issues!

With regards to the multi-component flash algorithms, we would absolutely be interested in incorporating a DE flash algorithm into Clapeyron.jl! A flash algorithm is really the only missing component to our package (aside from perhaps parameter estimation modules); we've been trying to find ways to implement either the RAND and HELD algorithms; the issue is finding the time to do these things... if you are willing and able, we would be extremely grateful if you shared your implement with us!

The speeds you mention are consistent with the values for a (very) basic RR implementation I've developed which is in no way robust... any improvement is more than welcome!

[denbigh]

Great. I don't have a lot of experience with pull requests, but I will try to educate myself. :) If I can't find the time in the next week, I could also send a Jupyter notebook with the code/tests.

[denbigh]

And yes, I can see HELD being a pain to implement (though in the long run, it is probably better than a DE implementation...)

[pw0908]

If you need any help navigating around GitHub / PR, feel free to dm us on Zulip!

[denbigh]

Sorry to bother about this again, but I'm still having some issues with SAFTgammaMie. In particular, mixtures of water and various species (e.g. methanol, CO2, ethanol) give NaN errors for most functions, especially at higher pressures (e.g. 100 bar) and when the water fraction is >85%. For example:

using Clapeyron

#Specify H2O Mole Fraction (error above ~x_water > 0.85)
x_water = 0.9

#Define Water-X mixture, where X may be various options
system = SAFTgammaMie(["water","carbon dioxide"])
#system = SAFTgammaMie(["water","methanol"])
#system = SAFTgammaMie(["water","ethanol"])

p = 100e5                                                    #Error for p > 10e5 Pa
T = 323.15                                                   #Error for T < 70C
x = [x_water,1-x_water]

println(gibbs_free_energy(system, p, T, x))
println(entropy(system, p, T, x))
println(fugacity_coefficient(system, p, T, x))
println(chemical_potential(system, p, T, x))

This prints various NaN's. If we replace SAFTgammaMie with other state functions (PR, PCSAFT) everything works well.
I'm particularly interested in looking at SAFTgammaMie predictions in pressurized CO2/water mixtures, so any comments/fixes would be much appreciated.
Thanks!

[pw0908]

This is just an example of how much tuning we have left to do... in short, the reason this is happening is because our initial guess for the volume is too small. The one we provide is generally ok for lower pressures. But, at higher pressures, we need to make the initial guess bigger. You can do this from the outside by just adding:

function Clapeyron.x0_volume_liquid(model::Clapeyron.SAFTgammaMieModel,T,z)
      v_lb = Clapeyron.lb_volume(model,z)
      return v_lb*2.0
end

I've tried this for your case and it works. N.B. Clapeyron.lb_volume(model,z) is a function that determine the volume at which the packing fraction of your system is 1. (which effectively leads to a bunch of log(0) problems). It's a good reference point to define your liquid volume.

This is also an example of a feature we hope to take more-advantage of: extensible methods. We hope, if users run into similar issues in the future, it will be easy enough for them to modify the underlying methods to suit their needs.

[longemen3000]

Weird, it seems that the volume solver is failing for the liquid phase, can you try using Clapeyron.volume_compress(model,p,T,x, V0=guess) (that's the internal volume solver, you can provide a guess until it converges)? , maybe our initial guessings for SAFT gamma Mie are not right and the value provided (Clapeyron.x0_volume(model,P,T,x,phase=:l)) is too far away from the solution

[longemen3000]

As a stopgap, you can do:

v = volume... #obtaining the volume from the specified pressure, by modifying the function or using the underlying solver
prop =VT_entropy(v, T, x) #calling the property with a VTn specification 

[denbigh]

Thanks for the responses. For your information, there are a few edge-cases where no initial guess seems to work with the solver. I've given an example below for water-methanol at 25C and 50 bar. I'm not sure why the solver fails for this case. I'm hesitant to guess a volume in any other way, as I suspect in condensed phases most properties are strongly dependent on the volume. However, for my own application the work-around of trying a few initial guesses until you get convergence seems to work most of the time. Thanks!

for multfactor = LinRange(1.25,2.5, 100)
    function Clapeyron.x0_volume_liquid(model::Clapeyron.SAFTgammaMieModel,T,z)
          v_lb = Clapeyron.lb_volume(model,z)
          return v_lb*multfactor
    end

    #Specify H2O Mole Fraction 
    x_water =0.85

    #Define Water-methanol mixture
    system = SAFTgammaMie(["water","methanol"])

    p = 50e5                           
    T = 298.15                        
    x = [x_water,1-x_water]

    println(Clapeyron.volume(system, p, T, x))
    
end

[pw0908]

After a bit of troubleshooting, I think I found the culprit. This is a general problem for most systems with a lot of association going on. Problems tend to occur when we exceed the maximum number of iterations required for the association fractions in SAFT-gamma Mie to converge (500). For my personal uses, I've raised that to 10000 (way over kill, but I'm dealing with systems which have a lot of association).

@longemen3000 , should we raise this limit in the default code? I've found increasing it to 1000 iterations is sufficient to deal with the above problem.

@denbigh For a quick fix, just paste this code at the beginning of your script (it increases the max iterations to 1000)::

function Clapeyron.X(model::Union{Clapeyron.SAFTModel,Clapeyron.CPAModel}, V, T, z,data = nothing)
    _1 = one(V+T+first(z))
    tol = model.absolutetolerance
    X_ = Clapeyron.PackedVectorsOfVectors.packed_ones(typeof(_1),length(@Clapeyron.sites(i)) for i ∈ @Clapeyron.comps)
    idxs = Clapeyron.indices(X_)    
    X0 = X_.v
    bv = model.params.bondvol.values
    nn = length(bv.values)
    if isone(nn)
        xia,xjb = Clapeyron.X_exact1(model,V,T,z,data)
        i,j = only(bv.outer_indices) 
        a,b = only(bv.inner_indices) 
        #in the case that i = j, a = b, this does assignment twice
        #we do xia last because xjb is calculated from xia
        X_[j][b] = xjb
        X_[i][a] = xia
        return X_
    end
    ρ = Clapeyron.N_A/V
    if data === nothing
        _Δ = @Clapeyron.f(Clapeyron.Δ)
    else
        _Δ = @Clapeyron.f(Clapeyron.Δ,data)
    end  
    fX(out,in) = Clapeyron.X!(out,in,idxs,_Δ,model.sites,ρ,z)
    Xsol = Clapeyron.Solvers.fixpoint(fX,X0,Clapeyron.Solvers.SSFixPoint(0.5),atol=tol,max_iters = 1000)
    return Clapeyron.PackedVofV(idxs,Xsol)
end

Sorry if it makes your code look ugly... this will hopefully only be a temporary fix.

[denbigh]

Fantastic, thanks! I just want to confirm that with both changes (change initial volume guess to 2*lb_volume, and increase iterations to 1000) everything works, and I get no NaN results for a broad sweep of conditions.

[longemen3000]

ok, i'm gonna bump the iteration limit, increase the x0_volume_liquid on SAFTgammaMie, and release a new version