Ternary LLE

[tlorance]

I am interested in using Clapeyron.jl to correlate experimental ternary LLE data. I can figure out how to construct a model (at least in some cases*), but I cannot figure out how to call "LLE" with three or more components: I don't think this is a bug, just a lack of documentation. Could someone please help?

*I did try to set up a NRTL model for water, acetone, and chloroform, but an error was thrown:
ERROR: Missing values exist ∈ single parameter Tc: Union{Missing, Float64}[647.13, 508.2, missing].
I assume that this means that the critical temperature of some component (chloroform, since it is last?) is missing, but I am suspicious I have misinterpreted this as the critical temperature of all of these compounds are well known. If it is simply missing from the code database, can I add it so that I can use NRTL?

[longemen3000]

It is a problem on how our current activity models are defined, try using NRTL(["water", "acetone", "chloroform"] puremodel = BasicIdeal)

[pw0908]

Just to add to this, the LLE function can only model LLE for a binary mixture. If your goal is to model LLE for a ternary mixture, I'd recommend using a flash algorithm instead. Here you do need to specify a non-ideal gas equation for the vapour phase. As you pointed out, we need to specify the critical point for all species. This can be done using the pure_userlocations optional argument. As an example, given chloroform is the one with the missing critical data here, we can define it in a csv file as such:

Clapeyron Database File,,,,,
Critical Single Parameters,,,,,
species,DIPPR Number,Tc,Pc,Vc,acentricfactor
chloroform,74828,537.,5.37,2.433e-4,0.218

I named this file example.csv. With this, you can construct your model using NRTL:

julia> model = NRTL(["water","ethanol","chloroform"];pure_userlocations=["example.csv"])
NRTL{PR{BasicIdeal, PRAlpha, NoTranslation, vdW1fRule}} with 3 components:
 "water"
 "ethanol"
 "chloroform"
Contains parameters: a, b, c, Mw

With this, you can now perform a flash calculation to obtain the LLE split:

julia> (x,n,G) = tp_flash(model,1e5,298.15,[0.4,0.2,0.4],RRTPFlash(equilibrium=:lle))
([0.8297322420540518 0.15921199667918248 0.011055761266765664; 0.03308014486558528 0.23482617035706146 0.7320936847773533], [0.3821558182414502 0.07332942819501094 0.005092032440156436; 0.0178441817585497 0.12667057180498906 0.39490796755984364], -13.173873832905187)

where x is the mole fraction in each phase, n is the molar amount in each phase and G is the gibbs free energy.

I hope this helps!

[tlorance]

This looks great. One question: where should I locate the "example.csv" file relative to my Julia installation? I am new to Julia, and I cannot find the Clapeyron.jl installation, so I assume its tucked away inside my Julia installation somewhere, but since I can't find it, it is not obvious to me where I should put files that Clapeyron needs to use during computation.

[pw0908]

If you start julia from your command line, the file just needs to be in the directory in which you started julia (similarly if you're using Jupyter notebook). If you start julia from the app, I think the file will have to be in your root directory.

[tlorance]

Thank you, I was able to get it to work using an absolute directory address for the example file. I am sorry to continue to be a bother, but I have a couple more questions:
First, is there any recourse when the computation fails to converge? For example, I set up the model as: model = NRTL(["water", "acetic acid", "chloroform"];pure_userlocations=["/Users/..."]) and the tp_flash command converged for many of our experimental bulk mixtures, but not for some: (x,n,G) = tp_flash(model,1e5,298.15,[0.417380784546791, 0.248249161677283, 0.334370053775926],RRTPFlash(equilibrium=:lle)) gave the starting mole fractions back for both phases and gave NaN back for n and G. I tried switching to MichelsenTPFlash, but it also failed, and I couldn't figure out how to even attempt DETPFlash as it apparently has a different signature from the other two. Would it help to use a different mixing rule? And, if so, how would I do that? I presume it would be an option in the model construction command...
Second, how would I go about verifying that the solution given by the flash algorithm is a stable and consistent solution. I see in the src directory that there are consistency and stability routines, but I can't figure out how to use them.
Thank you very much for all of your time.

[pw0908]

In this case, it means the default initial guesses we use didn't converge to the correct solution. In RRTPFlash() (and MichelsenTPFlash()), this can be done by specifying either K0 (initial partition coefficient). You can find more information on this in our documentation: https://clapeyronthermo.github.io/Clapeyron.jl/dev/properties/multi/#Clapeyron.RRTPFlash. As an example, I managed to converge the solution using:

julia> model = NRTL(["water","acetic acid","chloroform"];pure_userlocations=["example.csv"])
NRTL{PR{BasicIdeal, PRAlpha, NoTranslation, vdW1fRule}} with 3 components:
 "water"
 "acetic acid"
 "chloroform"
Contains parameters: a, b, c, Mw

julia> (x,n,G) = tp_flash(model,1e5,298.15,[0.417380784546791, 0.248249161677283, 0.334370053775926],RRTPFlash(equilibrium=:lle,K0=[1000.,10.,0.0001]))
([0.7375268155374085 0.24461033837441734 0.017862846088174216; 0.031865826087897405 0.25263097709477483 0.7155031968173279], [0.4029238173928658 0.13363491229775057 0.009758785692563492; 0.01445696715392478 0.11461424937953242 0.324611268083363], -13.054544125517769)

As for checking if the solution is consistent, you can just verify using the activity_coefficient:

julia> activity_coefficient(model,1e5,298.15,x[1,:]).*x[1,:]
3-element Vector{Float64}:
 0.7664098218637342
 0.21940560040722923
 0.7558330974181812

julia> activity_coefficient(model,1e5,298.15,x[2,:]).*x[2,:]
3-element Vector{Float64}:
 0.7664098218400759
 0.21940560040458024
 0.7558330973962782

As for stability, we do have the isstable function. Unfortunately, when I tested it in this case, I realised we developed it with full equations of state in mind. Technically the definition of stability in an activity coefficient model is simpler than this, but we haven't coded it up yet. Nevertheless, you can check if there are any nearby, stabler solutions using tpd:

julia> Clapeyron.tpd(model,1e5,298.15,x[1,:])
(Vector{Float64}[], Float64[], Symbol[], Symbol[])

julia> Clapeyron.tpd(model,1e5,298.15,x[2,:])
(Vector{Float64}[], Float64[], Symbol[], Symbol[])

So long as these solutions are empty, we can safely assume the phases are stable.

[tlorance]

I am getting the hang of it. However, I have a new question: where do you get the info for the compounds? I have no trouble finding Tc, Pc, Vc, etc., but the DIPPR numbers I can find (DIPPR IDs from the DIPPR Database) don't match the ones in the data block you gave me above. Also, it turns out you do have an entry for Trichloromethane, but its values don't match the data block you gave me for chloroform (they are similar, but different – including the DIPPR Number). Which values matter and which do not?
Thank you very much!

[pw0908]

For chloroform, I honestly just took the first result from NIST. They gave a range of values, but since LLE calculations do not depend on those calculations, I just picked a random one. The database we have from Clapeyron comes from a large database where they have averaged across a range of experimental values. As such, those are more reliable.

As for why the DIPPR numbers don't line up: I just copied the first row of our database file and only replaced the parameters. The DIPPR value here corresponds to methane. You don't actually have to include the DIPPR number within your custom parameters.

We initially had the DIPPR number there to avoid issues like this where one species can have two names. Thankfully, @longemen3000 is working on an update where this shouldn't be a problem anymore.

[tlorance]

OK, that makes sense. If I want to select a different mixing rule (such as HV, for example), how do I format that option? It seems, having looked into the code itself, that mixing rules are options for the EoS, but for NTRL the EoS is an option for NRTL. I did figure out how to change the EoS that NRTL uses, but not how to pass options to said EoS.

[longemen3000]

when using advanced mixing rules, the main EoS is the cubic. that is:

model = PR(["water","ethanol"],activity = NRTL, mixing = HVRule)

that will use the PR model, using the Huron-Vidal rule (That itself calls the NRTL model for obtaining gE).

When using cubics with advanced mixing rules (or any EoS in general defined in terms of helmoltz energy), the LLE problem is solved by solving the equality of chemical potentials (μi(p,T,x'') = μi(p,T,x')). that calculation may involve (or not) mixing rules

When using activity models, it always use the activity model for the liquid phase (or phases) and the pure model for saturation and properties of the gas phase. in LLE this is just γi(T,x'')*xi'' = γi(T,x')*xi',so no pure model is used when calculating LLE with just activity models (those are used when a gas phase appears)

you can also just prebuild the model before passing it as an option:

act = NRTL(["water","ethanol"])
model = PR(["water","ethanol"],activity = act , mixing = HVRule)

[tlorance]

So, just to be clear, model = NRTL(["water", "acetone", "dichloromethane"];puremodel=PR) would be equivalent to calling model = PR(["water", "acetone", "dichloromethane"], activity=NTRL, mixing=vdW1fRule) except that it would solve different equations to find the LLE? But if you want any mixing rule other than vdW1f, you must use the second syntax? Is there any practical difference in the two method of solving for LLE?

[pw0908]

There seems to be a bit of a mix up here. The two models you describe here are different. For context:

• Activity coefficient models can be used completely on their own (no cubic model needed) to model LLE behaviour. However, to model VLE behaviour through Raoult's law, you need another model to estimate the saturation pressure of each component. For simplicity, we decided to make the default here PR.
• In cubic equations of state, you need mixing rules to model mixtures. There's a class of mixing rules which pair up the cubic with an activity coefficient model called G^E mixing rules. The goal of these mixing rules is to set the excess gibbs free energy of both models to be equal to each other. While it may sound like this would mean that the cubic model should have the same predictions as the activity model for things like LLE, this is rarely the case.

To best visualise this, take a look at the Txy diagram in this example notebook: https://github.com/ClapeyronThermo/Clapeyron.jl/blob/master/examples/activity_models.ipynb. You can see that, while the end points of the VLE envelopes may be the same, the bubble and dew curves are not.

If you'd like to learn more on this, we have an introductory course: https://github.com/ClapeyronThermo/introduction-to-computational-thermodynamics. Part 2 is most-relevant here.

[pw0908]

As an aside, I've made some example notebooks for plotting ternary phase diagrams on a separate branch: https://github.com/ClapeyronThermo/Clapeyron.jl/blob/docs-update/examples/ternary_phase_diagrams.ipynb.

[tlorance]

Those posts were helpful. I am, unfortunately, stuck with an optimization problem. I am trying to use NRTL to model water-acetone-chloroform (which works just fine) and water-acetone-dichloromethane (which fails at every bulk mole fraction vector I try, no matter what K0 vector I start RRTP with). My models are: model = NRTL(["water","acetone","dichloromethane"];pure_userlocations=["etc"]) and model2=NRTL(["water","acetone","chloroform"];pure_userlocations=["etc"]) with files

Clapeyron Database File,,,,,
Critical Single Parameters,,,,,
species,DIPPR Number,Tc,Pc,Vc,acentricfactor
dichloromethane,1511,510.,6.08,1.93e-4,0.199

and

Clapeyron Database File,,,,,
Critical Single Parameters,,,,,
species,DIPPR Number,Tc,Pc,Vc,acentricfactor
chloroform,1521,536.4,5.472,2.39e-4,0.221902

The command I run is (x,n,G) = tp_flash(model,1e5,298.15,z0,RRTPFlash(equilibrium=:lle,K0=[10.,1000.,0.1])) with z0=[0.25,0.35,0.4] (although I tried it with many starting points and several sets of K-values). The Rachford-Rice method converges easily for model2 but will not converge for model. I tried giving it a starting x-vector (since I have experimental data which gives two phases), but it would not accept starting phase compositions; with x0=[0.4675,0.0895,0.443], it gives the error:

(x,n,G) = tp_flash(model,1e5,298.15,z0,RRTPFlash(equilibrium=:lle,x0=x0))
ERROR: invalid specification of initial points
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] MichelsenTPFlash(; equilibrium::Symbol, K0::Nothing, x0::Vector{Float64}, y0::Nothing, v0::Nothing, K_tol::Float64, ss_iters::Int64, nacc::Int64, second_order::Bool, noncondensables::Nothing, nonvolatiles::Nothing)
   @ Clapeyron ~/.julia/packages/Clapeyron/jVi2M/src/methods/property_solvers/multicomponent/tp_flash/Michelsentp_flash.jl:69
 [3] RRTPFlash(; equilibrium::Symbol, K0::Nothing, x0::Vector{Float64}, y0::Nothing, v0::Nothing, K_tol::Float64, max_iters::Int64, nacc::Int64, noncondensables::Nothing, nonvolatiles::Nothing)
   @ Clapeyron ~/.julia/packages/Clapeyron/jVi2M/src/methods/property_solvers/multicomponent/tp_flash/RachfordRicetp_flash.jl:59
 [4] top-level scope
   @ REPL[35]:1

[pw0908]

A few things to note:

• The reason model is failing is because Clapeyron doesn't have any binary parameters for dichloromethane. You can check this by doing (@longemen3000 we should probably throw an error for this?):

julia> model.params.a
3×3 PairParam{Float64}(["water", "acetone", "dichloromethane"]) with values:
 0.0    0.054  0.0
 6.398  0.0    0.0
 0.0    0.0    0.0

Since a 0 parameter implies that the two species have almost ideal interactions, you're not going to get a phase split. The way to resolve this would be to add the parameters manually just like you did for the cubic. In the case of activity coefficient models, the binary parameters are not symmetric. As such, if you provide parameter a_ij, you also need to provide a_ji. An example of how this file should look is shown below:

Clapeyron Database File,,,,,
NRTL Unlike Parameters,,,,,
species1,species2,a,b,c,source
chloroform,water,-7.352,3240.7,0.2,
water,chloroform,8.844,-1140.1,0.2,

There might be some NRTL parameters for your system online. Bear in mind that a and c are dimensionless, b is in kelvin. Different journals use different units (sometimes c is also called alpha).

• Small note: you can store both the chloroform and dichloromethane parameters in the same csv file.
• If you want to use x0 to define your initial conditions, you must also specify a y0 (or the composition of the second liquid phase). I personally prefer to use K0 since that's what the solver uses anyways.

[pw0908]

If you're stuck looking for parameters, there is also UNIFAC you could consider (all the groups seem to be present here) or SPTNRTL.jl which @longemen3000 is developing.

[tlorance]

I may need to use another method, since I can't seem to find coefficients for dichloromethane. I guess there is an opportunity for someone to do some research... There are a (very) few ternary systems with dichloromethane and one or the other of my desired components, so I could attempt to compute a_ij and a_ji, b_ij and b_ji, c_ij and c_ji from the data. If I were to do that, does "i" refer to the row in model.params.a (etc.) or does it refer to the columns?

[longemen3000]

looking at spt-NTRL and it's supplementary info (that is technically a normal NRTL, but the parameters were generated via a transformer neural network, paper: https://doi.org/10.1016/j.fluid.2023.113731),it seems that the training data includes those three molecules, so i suppose that the predictions for this specific system are good enough. SPTNRTL.jl is not registered yet, but i obtained the parameters from that specific subsystem:

sm = ["O", "CC(C)=O", "ClCCl"] #smiles for water, acetone, dicloromethane
#a0,a1,t0,t1,t2,t3 = SPTNRTL.spt_NRTL_params(sm) #TODO on my part, better api for this that don't require smiles
a0,a1,t0,t1,t2,t3 = ([0.0 0.228632 0.257017; 0.228632 0.0 0.244773; 0.257017 0.244773 0.0], [0.0 0.000136 0.000187; 0.000136 0.0 0.000193; 0.000187 0.000193 0.0], [0.0 13.374756 13.426971; 16.081848 0.0 18.17276; 13.804535 19.355072 0.0], [0.0 -415.527344 583.285156; -88.8125 0.0 -926.445312; 1163.712891 -579.439453 0.0], [0.0 -1.913689 -2.213547; -2.930901 0.0 -2.844528; -2.846214 -3.284027 0.0], [0.0 0.00153 0.004608; 0.005305 0.0 0.001674; 0.004722 0.004792 0.0])
model = aspenNRTL(["water", "acetone", "dichloromethane"],puremodel = BasicIdeal)
model.params.a0 .= a0
model.params.a1 .= a1
model.params.t0 .= t0
model.params.t1 .= t1
model.params.t2 .= t2
model.params.t3 .= t3

[tlorance]

Thank you. I was going to try (and would still like to try) regressing my data set to obtain the necessary values; I'm going to have to work my way up to that, though, as I don't immediately see what the appropriate estimator object would be.

[pw0908]

A few recommendations, if you'd like to do this in Clapeyron:

1. If you can, only use two-component experimental data (even VLE would be suitable for this). Trying to regress parameters with ternary LLE will be computationally intensive and may not result is the most generalised set of parameters.
2. If you need guidance on how to fit such parameters in Clapeyron, we have a notebook which shows how to do it in the case of activity coefficient models: https://github.com/ClapeyronThermo/Clapeyron.jl/blob/master/examples/parameter_estimation.ipynb.

If your work is time-sensitive, I do recommend you consider using either SPT-NRTL or UNIFAC. Parameter estimation can be time consuming on its own.

[tlorance]

Yes, I had already looked at that notebook and I follow it pretty well, having used similar systems in Python for other applications. I also intended to switch to fitting two-component systems (at least at first) so they would be more computationally tractable. However, for DCM-water, the available binary data is mutual solubility and I'm not sure how to set that up with respect to the estimator function. I suppose that the DCM-acetone binary would be easier since there I would be looking at bubble and dew points.

[tlorance]

New issue: I am attempting to use UNIFAC, but the group parameter for water seems to be missing. I tried model = UNIFAC(["water","acetone","dichloromethane"];pure_userlocations=["etc"]) and obtained the error ERROR: Missing values exist ∈ single parameter groups: Union{Missing, String}["[\"H2O\"=>1]", missing, missing].

[tlorance]

Also, using the spt-NTRL model and parameters you so kindly provided, I am getting the error ERROR: type BasicIdeal has no field components.

[pw0908]

A lot to answer here:

• I would recommend trying the Acetone-DCM system first if you already have the data, just to get a grip of parameter estimation. The computations should be quite fast. As for performing the water-DCM regression, since it is a binary, you can now switch back to the LLE function. You would just need to define a separate function for the exact compositions you wish to fit. This shouldnt be much more computationally intensive. I do recommend you avoid using data points close to the UCST.
• The error you're getting for UNIFAC is actually that the acetone and dichloromethane groups haven't been defined. You can define them as such:

julia> model = UNIFAC(["water",("acetone",["CH3CO"=>1,"CH3"=>1]),("dichloromethane",["CH2CL2"=>1])];pure_userlocations=["example.csv"])
UNIFAC{PR{BasicIdeal, PRAlpha, NoTranslation, vdW1fRule}} with 3 components:
 "water": "H2O" => 1
 "acetone": "CH3CO" => 1, "CH3" => 1
 "dichloromethane": "CH2CL2" => 1
Group Type: UNIFACDortmund
Contains parameters: A, B, C, R, Q

• The SPT-NRTL bug is mainly because the pure model was choosen to be BasicIdeal. You can fix this by just using PR again with the additional parameters you've defined.

[longemen3000]

The aspenNRTL bug is weird, what Clapeyron version are you using?

[tlorance]

@longemen3000 How do I find out? I have been unable to find the actual package location and I don't know how to request a version number in Julia. I installed it last week, so it should be pretty current.

[pw0908]

The bug is unrelated to your Clapeyron version. The fixed I mentioned previously should work?

If you want to check the Clapeyron version youre using, paste this command into Julia:

] status Clapeyron

[tlorance]

@pw0908 I am working through correlating my data with UNIFAC. It is going well, but I chanced upon a discovery that you might want to be aware of: the existing tests (as you noted for me above) for stability of the solution in an LLE are not foolproof. For example, if you use model = UNIFAC(["water",("acetone",["CH3CO"=>1,"CH3"=>1]),("dichloromethane",["CH2CL2"=>1])];pure_userlocations=["etc"]) and run

z1=[0.225,0.55,0.225]
K0=[8.74,0.26,0.03]
(x,n,G) = tp_flash(model,1e5,298.15,z1,RRTPFlash(equilibrium=:lle,K0=K0))

you will get the following test results:

activity_coefficient(model,1e5,298.15,x[1,:]).*x[1,:]
3-element Vector{Float64}:
 0.891532221995761
 0.5570887544557274
 0.2428284130952454

activity_coefficient(model,1e5,298.15,x[2,:]).*x[2,:]
3-element Vector{Float64}:
 0.8915322219984043
 0.5570887544411055
 0.24282841308151698

Clapeyron.tpd(model,1e5,298.15,x[1,:])
(Vector{Float64}[], Float64[], Symbol[], Symbol[])

Clapeyron.tpd(model,1e5,298.15,x[2,:])
(Vector{Float64}[], Float64[], Symbol[], Symbol[])

However, if you run

z1=[0.225,0.55,0.225]
K0=[1.536,0.9435,0.7848]
(x,n,G) = tp_flash(model,1e5,298.15,z1,RRTPFlash(equilibrium=:lle,K0=K0))

you get the following test results:

activity_coefficient(model,1e5,298.15,x[1,:]).*x[1,:]
3-element Vector{Float64}:
 1.0178459626398992
 0.509345369957719
 0.2771883900103174

activity_coefficient(model,1e5,298.15,x[2,:]).*x[2,:]
3-element Vector{Float64}:
 1.0191262487738564
 0.5091619606408663
 0.27687475194617833

Clapeyron.tpd(model,1e5,298.15,x[1,:])
(Vector{Float64}[], Float64[], Symbol[], Symbol[])

Clapeyron.tpd(model,1e5,298.15,x[2,:])
(Vector{Float64}[], Float64[], Symbol[], Symbol[])

The only way I could tell the second set was inconsistent was the sudden discontinuity in my phase diagram; there are an entire set of solutions similar to the second that exist in that region, but they are not continuous with the low-acetone tie-lines and they have higher (less negative) excess energies of mixing.

[pw0908]

This behaviour has more to do with the fact that the Rachford-Rice algorithm (RRTPFlash) is a local optimiser. It obtains the first phases close to your initial guess which satisfy chemical equilibrium (or as close as it can in your case). As such, if you try two different initial guesses, there's a chance it'll find two different solutions.

Algorithms like MichelsenTPFlash and DETPFlash are global optimisers which try to find the globally stable solutions. Unfortunately, DETPFlash is very slow and MichelsenTPFlash currently doesn't work for Activity coefficient models in LLE. There's a minor bug which we will try to fix soon.

My only recommendation for this case is to find one globally stable solution and then re-use this solution as the initial guess for the next iteration. You can see examples of this in the notebooks I shared for the ternary phase diagrams.

[tlorance]

Wow, you weren't kidding about DETPFlash. I ran it time-unrestricted and 24 hours later it still hadn't finished. I will try rerunning it with a time limit just to see what it does, but I guess I will be using RRTPFlash and human judgement.

[tlorance]

I find that I am unable to get near the plait point with RRTPFlash; the closer I get, the more it wants to diverge into alternate solutions and I cannot tell which, if any, are legitimate. I had been slightly varying the feed composition to improve convergence (since any feed composition along a tie line will give the same endpoints), but this alters the overall excess free energy of mixing so I can't tell if the solution is better, worse, or the same (by that criterion). The only criterion I have is consistency with previous results, and I am now worried that if I followed a suboptimal solution branch early on, I am just perpetuating that error by choosing solutions that are consistent with it. Is this a common problem with LLE, or is my system uniquely ill-conditioned?

[tlorance]

@pw0908 Am I doing this wrong: model = aspenNRTL(["water", "acetone", "dichloromethane"],puremodel = PR)? When I tried this, I still receive ERROR: type BasicIdeal has no field components.

[longemen3000]

I find that I am unable to get near the plait point with RRTPFlash; the closer I get, the more it wants to diverge into alternate solutions and I cannot tell which, if any, are legitimate. I had been slightly varying the feed composition to improve convergence (since any feed composition along a tie line will give the same endpoints), but this alters the overall excess free energy of mixing so I can't tell if the solution is better, worse, or the same (by that criterion). The only criterion I have is consistency with previous results, and I am now worried that if I followed a suboptimal solution branch early on, I am just perpetuating that error by choosing solutions that are consistent with it. Is this a common problem with LLE, or is my system uniquely ill-conditioned?

@pw0908 is likely a symptom of the bug we were talking earlier?

@pw0908 Am I doing this wrong: model = aspenNRTL(["water", "acetone", "dichloromethane"],puremodel = PR)? When I tried this, I still receive ERROR: type BasicIdeal has no field components.

It works on the master version of Clapeyron (so maybe something is wrong in the current version). i'm adding a test

Wow, you weren't kidding about DETPFlash. I ran it time-unrestricted and 24 hours later it still hadn't finished. I will try rerunning it with a time limit just to see what it does, but I guess I will be using RRTPFlash and human judgement.

This requires a fix on our part. i think you are the first on trying LLE with Differential evolution. on MichelsenTPFlash we cache a lot of things when activity models are involved. such caching is not done in DETPFlash, but that would be easy to solve.

I will try to address all those concerns ASAP, and release a new version. then i will ask you to update and try again with DETPFlash

in the meantime, i will add a testsuite just with the errors reported in this issue. feel free to write any code that you should expect to work but it doesn't

[tlorance]

It shouldn't matter, but I am working on a MacBook Pro; I installed Julia specifically to run Clapeyron, so the Julia installation should be up-to-date.

[pw0908]

In terms of actually tracing the LLE region using tp_flash, I don't know if you've taken the time to look at the example notebook but the best way to remain on the same stable branch as your initial point is to:

1. Use the current solution as the initial guess to the next iteration
2. Adapt your starting composition as you trace the two phase region such that it is always within the two phase region.

I went ahead an adapted my code for your system:

model = UNIFAC(["water",("dichloromethane",["CH2CL2"=>1]),("acetone",["CH3"=>1,"CH3CO"=>1])])

# Define initial conditions
T = 298.15
p = 1e5
z0 = [0.5,0.5,1e-10]

# Using the phase fraction to 'aim' for the plait point
ϕ = 0.75

# Give a decent initial guess
K0 = [1000.,0.0001,0.001]

# Define the tracing step size
dz = 1e-3

# initialise an empty vector which will store all of our points
x_LLE = zeros(500,6)

idxend = 500
for i in 1:500
     # Solve the flash
    (x,n,G) = tp_flash(model,p,T,z0,RRTPFlash(K0=K0,equilibrium=:lle))
    # Store the values
    x_LLE[i,1:3] = x[1,:]
    x_LLE[i,4:6] = x[2,:]
    
    # Use them as initial guesses in the next iteration
    K0 = [x[1,1]/x[2,1],x[1,2]/x[2,2],x[1,3]/x[2,3]]

    # Adapt your starting composition to be within the two-phase region
    z0[1:2] = x[1,1:2]*ϕ+x[2,1:2]*(1-ϕ)
    
    # Take a step further in the direction of the component that will close the two-phase region
    z0[3] = z0[3] + dz

    if x[1,1]-x[2,1] < 1e-3
        # If the two-phase region closes, give the index at which it closes and stop the loop
        idxend=i-1
        break
    end
end

# Store all the values into a single vector for plotting
x_LLE = vcat(x_LLE[1:idxend,1:3],reverse(x_LLE[1:idxend,4:6],dims=1));

The figure I got looks correct:

[tlorance]

The figure looks great; if I am reading the code correctly, tie lines would connect x[n,1:3] and x[n,4:6]. This is essentially what I tried to do, including predicting the next K0 from the previous round, stepping the acetone gradually, and so forth. The one thing I didn't do (and it seems to have made all the difference) is use the method you used to predict the next z0 values. Also, I didn't go down to 0.001-sized steps until I ran into trouble.
One question: what modification would you make to identify the plait point? Is it just the coordinates at index idxend?
Also, out of curiosity, is the plot from the Clapeyron routine (I haven't prioritized reading that part since I know how to plot ternary phase diagrams in R and python)?
Thank you very much!

[longemen3000]

Also, out of curiosity, is the plot from the Clapeyron routine (I haven't prioritized reading that part since I know how to plot ternary phase diagrams in R and python)?

There aren't any integrated plotting functions in Clapeyron (yet) but you call PyPlot directly in Julia (using PyPlot)

[tlorance]

So I copied your script into a file and ran it in Julia (using Clapeyron) and it is throwing two warnings and an error:

┌ Warning: Assignment to `K0` in soft scope is ambiguous because a global variable by the same name exists: `K0` will be treated as a new local. Disambiguate by using `local K0` to suppress this warning or `global K0` to assign to the existing global variable.
└ @ ~/Documents/Local_Projects/SURP 2019/Ternary Phase/Edwin Data/WaterDCMAcetone_UNIFAC.txt:29
┌ Warning: Assignment to `idxend` in soft scope is ambiguous because a global variable by the same name exists: `idxend` will be treated as a new local. Disambiguate by using `local idxend` to suppress this warning or `global idxend` to assign to the existing global variable.
└ @ ~/Documents/Local_Projects/SURP 2019/Ternary Phase/Edwin Data/WaterDCMAcetone_UNIFAC.txt:39
ERROR: LoadError: UndefVarError: `K0` not defined
Stacktrace:
 [1] top-level scope
   @ ~/Documents/Local_Projects/SURP 2019/Ternary Phase/Edwin Data/WaterDCMAcetone_UNIFAC.txt:23
 [2] include(fname::String)
   @ Base.MainInclude ./client.jl:478
 [3] top-level scope
   @ REPL[2]:1
in expression starting at /Users/ted/Documents/Local_Projects/SURP 2019/Ternary Phase/Edwin Data/WaterDCMAcetone_UNIFAC.txt:21

The warnings I'm not too worried about, since I think that treating the variables as local is correct. However, I can't figure out why it says that K0 isn't defined, especially since (if I am counting code lines correctly) it is claiming it isn't defined the third time K0 appears! I'm not sure what to do...

[longemen3000]

try wrapping all that code in a function:

function lle_solve()
model = UNIFAC(["water",("dichloromethane",["CH2CL2"=>1]),("acetone",["CH3"=>1,"CH3CO"=>1])])

# Define initial conditions
T = 298.15
p = 1e5
z0 = [0.5,0.5,1e-10]

# Using the phase fraction to 'aim' for the plait point
ϕ = 0.75

# Give a decent initial guess
K0 = [1000.,0.0001,0.001]

# Define the tracing step size
dz = 1e-3

# initialise an empty vector which will store all of our points
x_LLE = zeros(500,6)

idxend = 500
for i in 1:500
     # Solve the flash
    (x,n,G) = tp_flash(model,p,T,z0,RRTPFlash(K0=K0,equilibrium=:lle))
    # Store the values
    x_LLE[i,1:3] = x[1,:]
    x_LLE[i,4:6] = x[2,:]
    
    # Use them as initial guesses in the next iteration
    K0 = [x[1,1]/x[2,1],x[1,2]/x[2,2],x[1,3]/x[2,3]]

    # Adapt your starting composition to be within the two-phase region
    z0[1:2] = x[1,1:2]*ϕ+x[2,1:2]*(1-ϕ)
    
    # Take a step further in the direction of the component that will close the two-phase region
    z0[3] = z0[3] + dz

    if x[1,1]-x[2,1] < 1e-3
        # If the two-phase region closes, give the index at which it closes and stop the loop
        idxend=i-1
        break
    end
end

# Store all the values into a single vector for plotting
x_LLE = vcat(x_LLE[1:idxend,1:3],reverse(x_LLE[1:idxend,4:6],dims=1));
return x_LLE
end

that error is basically julia telling us that we really shouldn't use globals for everything, you can read more of that here: https://docs.julialang.org/en/v1/manual/variables-and-scoping/#on-soft-scope

[pw0908]

Answering your questions:

1. Identifying the Plait point is tough...if you want it approximately you can just average the two points at idxend. To improve the precision, you can reduce the step size and increase the number of points. However, the closer you get to the plait point, the more unstable the solvers will get and you are likely to run into errors (just the nature of critical points). Alternatively, you can use a numerical solver to solve for the conditions of criticality + pressure constraint to solve for the compositions of the plait point. This might be a bit challenging as you'd have to look into our source code. If it works, it would get you the exact location of the plait point.
2. On the K0 error, that's my bad. I coded this in Jupyter Notebook and didn't account for the for loop restrictions in julia (for loops are treated as a separate computational environment in julia; anything defined outside the for loop will not be defined inside). @longemen3000 correct me if I'm wrong. Either way, the fix is simply (add global in the for loop):

model = UNIFAC(["water",("dichloromethane",["CH2CL2"=>1]),("acetone",["CH3"=>1,"CH3CO"=>1])])

N = 500
T = 298.15
p = 1e5
z0 = [0.5,0.5,1e-10]

K0 = [1000.,0.0001,0.001]
x_LLE = zeros(N,6)
idxend = N
(x,n,G) = tp_flash(model,p,T,z0,RRTPFlash(K0=K0,equilibrium=:lle))
x_LLE[1,1:3] = x[1,:]
x_LLE[1,4:6] = x[2,:]

K0 = x[1,:]./x[2,:]
z0[1:2] = x[1,1:2]/2+x[2,1:2]/2
z0[3] += 2/N
for i in 2:N
    global z0, K0, x_LLE, idxend
    (x,n,G) = tp_flash(model,p,T,z0,RRTPFlash(K0=K0,equilibrium=:lle))
    x_LLE[i,1:3] = x[1,:]
    x_LLE[i,4:6] = x[2,:]
    K0 = x[1,:]./x[2,:]

    z0 = (x_LLE[i,1:3]+x_LLE[i,4:6])-(x_LLE[i-1,1:3]+x_LLE[i-1,4:6])/2
    if abs(x[1,1]-x[2,1]) < 1e-3
        idxend=i-1
        break
    end
end

x_LLE = vcat(x_LLE[1:idxend,1:3],reverse(x_LLE[1:idxend,4:6],dims=1));

[pw0908]

@tlorance sorry about that, we should've probably communicated before answering at the same time 😆. Either solution is fine.

[tlorance]

What is the format for supplementing UNIQUAC parameters? I have found parameters for dichloromethane, but I can't quite figure out the syntax required to supply them to the model.
Thank you very much.

[tlorance]

I have the same question about COSMO-SAC, although I'm having a harder time finding those parameters (and, in that case, acetone is also missing).

[pw0908]

As a general tip, the parameters for each equation of state is described in the docs; if you wish to input custom parameters, you can either feed them in through CSVs using the naming conventions described in the docs or directly within the code. Nevertheless, for your specific case:

• For UNIQUAC, you just need to do something similar to what I described earlier with NRTL, although you now need both like and unlike parameters:

For the unlike parameters (unlike.csv):

Clapeyron Database File,,,,,
UNIQUAC Unlike Parameters,,,,,
species1,SMILES1,species2,SMILES2,a,source
methanol,,benzene,,-123.69,
benzene,,methanol,,935.05,

For the like parameters (like.csv):

Clapeyron Database File,,,,
UNIQUAC Like Parameters,,,,
species,DIPPR Number,r,q,q_p,source
methanol,,1.43,1.43,1.43,
benzene,,3.19,2.4,2.4,

q and q_p should be the same unless you are using a special variant of UNIQUAC. With this, you can build the model like you did before (although now, you specify userlocations rather than pureuserlocations):

model = UNIQUAC(["methanol","benzene"]; userlocations=["like.csv","unlike.csv"])

As an aside, why are you looking at UNIQUAC? UNIFAC already has the parameters for dichloromethane. The approaches will be identical for species represented as a single group.

• COSMO-SAC is a bit trickier. Here you need sigma profiles as the parameters, but the type of profile you need depends on which version of COSMO-SAC you use. COSMO-SAC 02 only requires the full sigma profile, COSMO-SAC 10 requires three different profiles and COSMO-SAC dsp requires the same three profiles plus some pure parameters. In the case you need to use COSMO-SAC dsp, all the parameters can be found in this repository. It's unfortunately an uneasy repository to navigate through. Once you've found the relevant profile, you'll need to reformat it into a CSV compatible with Clapeyron (see the COSMOSAC10_like.csv and COSMOSACdsp_like.csv files for reference on how to format the CSVs). Once you've done all of that, just do as you were doing previously for all the other equations of state.

This is a bit convoluted but there is a copyright on those parameters meaning that we cannot store them locally in Clapeyron (the files are also massive, which would unnecessarily increase the size of Clapeyron). We do have some ideas on how to streamline the process using special parsers, but these are not high priority right now.

Finally, as this thread has now evolved beyond a single issue and the fact that it is now getting very long, would you be open to switching to email? Our addresses can be found in the docs (Pierre and Andrés).

[tlorance]

I thought that UNIFAC used group-contribution methods to estimate the UNIQUAC parameters, so I had assumed there was probably a difference between the UNIFAC-estimated parameters and the fitted parameters used in base UNIQUAC.
To be honest, I'm not as concerned with the precise model: I am trying to publish some sets of experimental ternary phase data with a quantum chemical solvation interpretation. However, having looked at the literature, I have observed that one does not publish ternary LLEs without fitting them to at least one, and preferably 2+, activity models.
I also thought that, since the environmental hazard of DCM is an up-and-coming topic, it would be helpful to have more activity coefficient data for it, so if it turned out that our data could contribute, all the better.
I could absolutely switch to email; I was thinking about that very problem with this thread myself. Would the two of you prefer I address you both when I write? I will do so the first time; then you will both have my email at least.

[pw0908]

• UNIFAC is was you might refer to as a heterogeneous group-contribution approach where, rather than obtaining one set of parameters to represent a species, it uses groups to represent different functional groups of a species (in a large alcohol, it doesn't make much sense to use the same parameters for the alkyl and hydroxyl groups?). UNIFAC reduces to UNIQUAC when species are represented by single groups. There might be a difference here as acetone is represented using two groups in UNIFAC.
• In terms of what is done in the literature, I believe I usually see UNIFAC, NRTL and Wilson as the usual suspects that people fit. Adding COSMO-SAC to the mix might be interesting since it is a purely predictive approach. You could also include other models such as predictive cubics and SAFT equations. At that point, it would become more of a comparative study rather than just presenting new data. This might supersede the need to fit new parameters?
• Feel free to address us both. This'll just speed things up on your end since youll only need to wait for one of us to reply.

[pw0908]

Conversation has moved to email so I'll close the issue