Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Something new that's not in the database #193

Open
1578jason opened this issue Aug 1, 2023 · 16 comments
Open

Something new that's not in the database #193

1578jason opened this issue Aug 1, 2023 · 16 comments

Comments

@1578jason
Copy link

For some complex compounds, many models do not seem to work, such as the SAFTgammaMie model for the calculation of n-tetracosane, how to construct the model and give the corresponding model parameters, can you give an example?

@pw0908
Copy link
Member

pw0908 commented Aug 1, 2023

Hi!
You can define new species in SAFT-gamma Mie using:

julia> using Clapeyron

julia> model = SAFTgammaMie([("tetracosane",["CH3"=>2,"CH2"=>22])])
SAFTgammaMie{BasicIdeal, SAFTVRMie{BasicIdeal}} with 1 component:
 "tetracosane": "CH3" => 2, "CH2" => 22
Group Type: unknown
Contains parameters: segment, shapefactor, lambda_a, lambda_r, sigma, epsilon, epsilon_assoc, bondvol

julia> crit_pure(model)
(857.6313015851466, 1.0961825768823908e6, 0.0022801353592358043)

Hope this helps!

@rwer78
Copy link

rwer78 commented Aug 8, 2023

Hi! I have a similar question about using the SAFT equation of state to define a component .
I want to do some calculations using SAFT-type equations of state, where the substances and the corresponding parameters are self-defined. There is no problem when I use PCSAFT model to define some non-associated substances, but there is an error when I customize associated compounds, can you help me correct and give correct examples? In addition, I also made a similar error when using SAFTVRMie(using our own params). Can you give me a SAFTVRMie example(using our own params) by the way?
mmexport1691518891382

@longemen3000
Copy link
Member

of course, for that, you need to specify:

  • the amount and type of sites
  • the assoc params in a recognized form

in this particular example,

model = PCSAFT(["methanol"];userlocations =(;
    Mw = [32.042],
    segment = [1.5255],
    sigma = [3.23],
    epsilon = [188.9],
    n_H = [1],
    n_e = [1],
    epsilon_assoc = Dict((("methanol","e"),("methanol","H")) => 2899.5),
    bondvol = Dict((("methanol","e"),("methanol","H")) => 0.035176))
    )

note how we specify the amount of sites (e is specified by n_e and H is specified by n_H) and how the association parameter is specified by a dict, containing a relation between a pair of component-sites and the value.
on SAFTVRMie:

model = SAFTVRMie(["methanol"];userlocations =(;
    Mw = [32.042],
    segment = [1.67034],
    sigma = [3.2462],
    epsilon = [307.69],
    lambda_a = [6.0],
    lambda_r = [19.235],
    n_H = [1],
    n_e = [1],
    epsilon_assoc = Dict((("methanol","e"),("methanol","H")) => 2062.1),
    bondvol = Dict((("methanol","e"),("methanol","H")) => 1.0657e-28))
    )

@rwer78
Copy link

rwer78 commented Aug 9, 2023

Thank you very much! I have a new question about it. Suppose I haven't found the PCSAFT model parameters for a new associated component, how do I build the model? Or do I need to use experimental data for regression, are there any examples of that?

@pw0908
Copy link
Member

pw0908 commented Aug 9, 2023

You have a few options. The best would be to build the model yourself using parameter estimation with experimental data (which we provide the tools for). If your species is lacking experimental data, you could consider group contribution approaches (such as gcPCSAFT which is also available in Clapeyron) where you can assemble species from groups.

What species are you trying to model?

@pw0908
Copy link
Member

pw0908 commented Aug 10, 2023

As for examples, you can check out this jupyter notebook: https://github.com/ClapeyronThermo/Clapeyron.jl/blob/master/examples/parameter_estimation.ipynb

@pw0908 pw0908 closed this as completed Sep 8, 2023
@1578jason
Copy link
Author

Hi!
I have some questions about the parameter estimation of binary system. I noticed that the parameter estimation in the notebook is to directly regress the final result of a model parameter. Can I calculate BIP by regression?

@pw0908
Copy link
Member

pw0908 commented Sep 20, 2023

Hi! Yes, you can indeed fit BIP using our regression tool. You'd have to fit the epsilon_12 (in SAFT equations) or the a_12 (in cubics) to get BIP (you can then convert this to a kij parameter). I realised there isn't a good example of this in the notebooks. In short, you can follow the example from the activity coefficient models in here: https://github.com/ClapeyronThermo/Clapeyron.jl/blob/master/examples/parameter_estimation.ipynb. Your parameters are now symmetric so you don't need to worry about the asymmetric keyword.

Im currently working on updating the notebooks and docs. Im going to include more relevant examples in these.

@1578jason
Copy link
Author

Hi! I'm sorry to trouble you. I noticed that when calculating gas-liquid equilibrium properties (such as bubble point function), T,x and P,x are usually specified to calculate the other two properties. Is it possible to specify the two variables of temperature and pressure to calculate the corresponding gas phase composition and liquid phase composition in gas-liquid equilibrium in Clapeyron?

@pw0908
Copy link
Member

pw0908 commented Oct 31, 2023

Hi! What you are describing is a flash algorithm. We support these in Clapeyron. Here is an example that will only look for VLE (the default method in Clapeyron makes no assumptions about the phase equilibrium type and is therefore quite slow):

julia> model = PCSAFT(["water","carbon dioxide"])
PCSAFT{BasicIdeal, Float64} with 2 components:
 "water"
 "carbon dioxide"
Contains parameters: Mw, segment, sigma, epsilon, epsilon_assoc, bondvol

julia> (x,n,G) = tp_flash(model,1e5,298.15,[0.5,0.5],MichelsenTPFlash(equilibrium=:vle))
([0.9995531870252333 0.00044681297476665416; 0.033441597210702265 0.9665584027892977], [0.48270815023652947 0.00021577667636994686; 0.017291849763470462 0.4997842233236301], -6.574034182311726)

julia> x
2×2 Matrix{Float64}:
 0.999553   0.000446813
 0.0334416  0.966558

x contains the composition in each phase.

@1578jason
Copy link
Author

Hi! I'm sorry to bother you again. I'm trying to do some phase equilibrium calculations with Clapeyron. According to the p-xy example in the notebook, I used PR EOS to calculate the bubble point pressure and gas composition of methane and ethane system at 280K (experimental data from the literature), and the results are shown below. All 15 data points have obtained the calculated results. At the same time, when I calculated directly using the bubble_pressure function, I got a different result (some points have inconsistent results and some points could not be calculated). Could you help me to check it why?

using Clapeyron
model=PR(["methane","ethane"])
T=280
xexp=[0.0081,0.0231,0.0335,0.0693,0.0996,0.1292,0.1528,0.1798,0.2044,0.2333,0.2668,0.2774,0.2913,0.3046,0.3109]
yexp=[0.0279,0.0744,0.1036,0.1864,0.2404,0.2818,0.3091,0.3336,0.3512,0.3654,0.3719,0.3711,0.3675,0.3609,0.3557]
pexp=[28.88,30.75,32.04,36.56,40.38,44.16,47.12,50.41,53.30,56.46,59.84,60.87,61.91,62.66,62.92]
X = Clapeyron.FractionVector.(xexp)

y = []
p = []
bub = []
v0 =[]
for j ∈ 1:15
    if j==1
        append!(bub, [bubble_pressure(model,280,X[j])])
        v0 = [log10(bub[j][2]),log10(bub[j][3]),bub[j][4][1],bub[j][4][2]]
    else
        append!(bub, [bubble_pressure(model,280,X[j];v0=v0)])
        v0 = [log10(bub[j][2]),log10(bub[j][3]),bub[j][4][1],bub[j][4][2]]
    end
end
append!(y,[append!([bub[i][4][1] for i ∈1:15])])
append!(p,[append!([bub[i][1] for i ∈ 1:15])])

println(y[1])
println(p[1])#[2.9090265171896922e6, 3.0874681588960933e6, 3.211008647049053e6, 3.634571295771231e6, 3.9899441569014466e6, 4.3328824450574415e6, 4.602135150580931e6, 4.90402726725854e6, 5.171510000245315e6, 5.47320087747928e6, 5.79869438646223e6, 5.89425338938065e6, 6.01238870391992e6, 6.115980641807861e6, 6.161079995839753e6]

bub1=bubble_pressure.(model,T,X[j] for j ∈ 1:15)

for i ∈ 1:15
    println(bub1[i][1])
end
#[2.9090265171896922e6, 3.0874681550907055e6, 3.211008643618842e6, 3.634571307084665e6, 3.9899441568940105e6, 4.332882436602619e6, 4.602135150581037e6, 4.904027267026698e6, 5.171510000248986e6, 5.473200877667004e6, NaN, NaN, NaN, 5.336779399714725e6, NaN]

@pw0908
Copy link
Member

pw0908 commented Nov 2, 2023

Hi Jason,

The reason for this because your mixture has a supercritical component (methane). Within Clapeyron, to generate initial guesses, we assume that we have an ideal mixture that obeys Raoult's law. Unfortunately, that also means that assume all components are subcritical. When that isn't the case, we make some basic approximations but it doesn't work very well frame from a species' critical point. At 280K, you are very far from methane's critical point (190K). This means our solvers will be using quite poor initial guesses.

The workaround for this is to trace the vle region by re-using the previous solution as the initial guess for the next iteration. These are much better guesses and are more likely to run into the correct solution.

Are you trying to run parameter estimation with these mixtures?

@pw0908 pw0908 reopened this Nov 2, 2023
@1578jason
Copy link
Author

Yes, thank you very much for your response and I believe Clapeyron is a very good tool to help me achieve this.

@pw0908
Copy link
Member

pw0908 commented Nov 3, 2023

If you are planning to do parameter estimation, I recommend you use tp_flash instead

@1578jason
Copy link
Author

Hi! I'm so sorry to bother you again. I have several questions when using flash algorithm to do calculations:

  1. The calculation seems not to converge to the correct result. Both n and G give NaN, which I guess is the reason for the initial value, but I don't know how to set a better initial value.
  2. When doing flash calculation, the designation of n also seems to have an impact on the calculation result, how to ensure that the composition of each calculation is in the two-phase region (VLE).
  3. Can you give me an example of visualizing VLE phase diagram with flash algorithm? (example of actually tracking VLE phase diagram,x-y phase diagram ). I have tried several times without success,because I don't know much about flash algorithm.
    Thank you very much.
model=PR(["benzene","heptane"])
(x,n,G) = tp_flash(model,0.2367e5,328.15,[0.1,0.9],MichelsenTPFlash(equilibrium=:vle))

@pw0908
Copy link
Member

pw0908 commented Nov 10, 2023

Hi, I must apologies for the delay in my reply, I have just come back from a conference.

I am a little bit confused: what is your objective here? Is it to just trace the phase envelope or to do parameter estimation.

If your goal is to trace the phase envelope, then you should use the bubble_pressure and dew_pressure as these functions will follow the phase envelope exactly. Flash algorithms can be used but you will need to follow the phase envelope manually.

If your goal is to do parameter estimation, then you'll want to fit to the composition in each phase at a given pressure and temperature. Here, the bubble and dew_pressure solvers are useless. This is why using flash algorithms is much better.

From what you're saying, you seem to be more interested in just tracing the phase boundary. In that case, please take a look at the example notebooks.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

4 participants