Phase diagrams?

[infinity77]

Hi,

thank you for making your code available. I am no expert and maybe I am asking very stupid questions, so please bear with me.

I deal with multi-components gas mixtures (with generally 8-10 components) and I would like to generate phase diagrams (pressure and temperature on x and y axes) for the mixtures. Do you know how I can go about leveraging your work to achieve that?

thank you in advance.

[gustavochm]

Hey!

You can generate phase diagrams (P-T projection) by using the flash function recursively. This function will attempt to solve the phase equilibria at given global composition, pressure and temperature and it requires initial guesses of the phases. Those can be generated by minimizing the tangent plane distance.

Hope this help!

[pferna6]

Hello,

I am also attempting to create phase diagrams and identify the critical points of binary mixtures using flash calculations but I am having difficulty with the code converging to a solution. Do you have any recommendations on how I can go about fixing this error?

[gustavochm]

Hi!

Can you explain further with type or problems are you having? Is the flash routine directly diverging (e.g., getting nan) or is it converging to a single phase point (beta = 0. or 1.)?

Critical points calculations are difficult to compute in the sense they usually require a really a good initial guess. If you are constructing the whole phase diagram, a simple hack would be to have more points near the critical zone and reuse the previous solution as the initial guess of the next point calculation. If you are having problem with initial guesses I would suggest using the tpd functions to get them.

Those are some general recommendations but I would need extra information to help you more.

Gustavo

[pferna6]

Hi Gustavo,

The flash routine was not diverging. As you mentioned, identifying the critical points (cp) can be quite difficult. Rather than identifying the cp directly. I am aiming to identify the isothermal phase diagrams for my system of interest to help approximate the parameters for the cp. Please correct me if I am mistaken but the flash calculations are for when the mixture has two phases (liquid vapor). This would be within the phase diagrams and would not help me identify the phase diagram curves. Do you recommend that I simply create a loop that will plot out the Pressure versus composition isothermal phase diagrams using Phasepy's bubble and dew functions?

[gustavochm]

Hey!

Yes you are right, the two phase flash calculation are mean to compute two phases in equilibria at a given global composition (z), temperature and pressure. The flash can be used for vapour-liquid and liquid-liquid equilibria.

If you are not getting nan and the flash calculation is returning just a single phase, this could be due to the global composition (if the global composition is not between the two phases, the flash will return a single phase).

If you are interested on the phase diagram at a single temperature it makes sense to compute it using bubble points calculations. As you don't know yet where the critical point is, I would suggest iterating from compositions from 0 to 1 and checking if the compositions of the phases are equal.

This is an example from a mixture of methane + decade which using SAFT-VR-Mie (sgtpy):

methane = component('methane', ms = 1.0, sigma = 3.752 , eps = 170.75,
                    lambda_r = 16.39, lambda_a = 6.)

decane = component('decane', ms = 3.0, sigma = 4.584 , eps = 415.19,
                    lambda_r = 20.92, lambda_a = 6.)

kij = -0.07948919865560719
Kij = np.array([[0, kij], [kij, 0]])
mix = mixture(methane, decane)
mix.kij_saft(Kij)
eos = saftvrmie(mix)

T = 71 + 273.15 #K

n = 500
# xf here is the final composition which was close to the critical point
xf = 0.85
x1 = np.linspace(0, xf, n)
X = np.array([x1, 1-x1])
T = 71 + 273.15 #K

Y = np.zeros_like(X)
P = np.zeros_like(x1)
vl = np.zeros_like(x1)
vv = np.zeros_like(x1)

i = 0 
P0 = 2.6e3
sol = bubblePy(X[:,i], P0, X[:,i], T, eos, full_output= True)
Y[:,i], P[i] = sol.Y, sol.P
# saving the volumes here is just to speed-up SAFT-VR-Mie in following iterations
vl[i], vv[i] = sol.v1, sol.v2

for i in range(1,n):
    sol = bubblePy(Y[:,i-1], P[i-1], X[:,i], T, eos, v0 = [vl[i-1], vv[i-1]], 
                   good_initial = True, full_output= True)
    Y[:,i], P[i] = sol.Y, sol.P
    vl[i], vv[i] = sol.v1, sol.v2

I hope that helps,

Gustavo

[Danlira1991]

Hello Gustavo, I'm trying to plot phase diagram of 9 components (c1,c2,c3,n-c4,ic4,n-c5,i-c5,c6,c7+) using bubblePy and dewPy. First I generete this 9 comp. pahse diagram in Hysys in order to compare with PhasePy. I try to compute point by point in bubblePy rotine but end up getting trivial solution, I mean my pressure answer is the pressure guessed. Then i try to use this 'for' loop in order to compute all point at once but get NaN solutions. Do you have any thoghts about it? How can I get a bubbleP and dewP line for this system using a 'for' loop?

[gustavochm]

Hey Danlira,

A trivial solution is usually obtained when the initial guess is poor (either in composition or pressure). I guess that for a 9 component mixture to generate a good initial guess is a challenging but I would still try to use the tpd functions to see if at least you can improve the composition initial guess.

In the case of the for loop, this needs for your first point to be computed successfully because it reuses the previous solution as the initial guess for the next iteration. If that's the case the good_initial=True is suitable, oversize I would recommend to change it back to false.

By the way which model/equation of state are you using? Also, if you don't mind sharing I can give a look to your files 😋

Gustavo

[Danlira1991]

Hey Danlira,

A trivial solution is usually obtained when the initial guess is poor (either in composition or pressure). I guess that for a 9 component mixture to generate a good initial guess is a challenging but I would still try to use the tpd functions to see if at least you can improve the composition initial guess.

In the case of the for loop, this needs for your first point to be computed successfully because it reuses the previous solution as the initial guess for the next iteration. If that's the case the good_initial=True is suitable, oversize I would recommend to change it back to false.

By the way which model/equation of state are you using? Also, if you don't mind sharing I can give a look to your files 😋

Gustavo

Hello Gustavo, Thank you for the anwser.

That is my file:

import numpy as np import matplotlib.pyplot as plt from phasepy import mixture, component, preos, cubic from phasepy.equilibrium import bubblePy, dewPx, bubbleTy

methane = component(name='methane', Tc=190.7, Pc=46.41, Zc=0.286, w=0.0115) ethane = component(name='ethane', Tc=305.43, Pc=48.84, Zc=0.279, w=0.0986) propane = component(name='propane', Tc=369.9, Pc=42.57, Zc=0.276, w=0.1524) i_butane = component(name='i_butane', Tc=408.05, Pc=36.48, Zc=0.282, w=0.1848) n_butane = component(name='n_butane', Tc=425.15, Pc=37.97, Zc=0.274, w=0.2010) i_pentane = component(name='i_pentane', Tc=460.35, Pc=33.34, w=0.2222) n_pentane = component(name='n_pentane', Tc=469.65, Pc=33.75, Zc=0.27, w=0.2539) n_Hexane = component(name='n_Hexane', Tc=507.85, Pc=30.32, Zc=0.266, w=0.3007) n7_plus = component(name='n7_plus', Tc=712.85, Pc=15.89, w=0.6530)

mix = methane + ethane+ propane +i_butane + n_butane+ i_pentane+n_pentane+n_Hexane+n7_plus Kij=np.array([[0, 0.00224137306213379, 0.00682878494262695, 0.0131134390830994, 0.0123046040534973, 0.017627477645874, 0.0179253816604614 ,0.0234740376472473, 0.0576416254043579], [0.00224137306213379, 0, 0.00125795602798462, 0.00457346439361572, 0.0040963888168335, 0.00741332769393921, 0.00760942697525024, 0.0114138126373291, 0.0382192730903625], [0.00682878494262695, 0.00125795602798462, 0, 0.00104051828384399, 0.000818967819213867, 0.00258338451385498, 0.00270050764083862, 0.0051419734954834 ,0.0260568857192993], [0.0131134390830994, 0.00457346439361572, 0.00104051828384399, 0, 0.0000133514404296875 , 0.000346183776855469, 0.000390052795410156, 0.00156527757644653, 0.0169016718864441], [0.0123046040534973, 0.0040963888168335, 0.000818967819213867, 0.0000133514404296875, 0, 0.00049513578414917, 0.000547230243682861, 0.00186634063720703, 0.0178416967391968],
[0.017627477645874, 0.00741332769393921, 0.00258338451385498, 0.000346183776855469, 0.00049513578414917,0, 0.0000012516975402832, 0.000439941883087158, 0.0124849081039429], [0.0179253816604614, 0.00760942697525024, 0.00270050764083862, 0.000390052795410156, 0.000547230243682861, 0.0000012516975402832, 0, 0.000393450260162354, 0.0122348666191101], [0.0234740376472473, 0.0114138126373291, 0.0051419734954834, 0.00156527757644653, 0.00186634063720703, 0.000439941883087158, 0.000393450260162354, 0,0.00828582048416138], [0.0576416254043579, 0.0382192730903625, 0.0260568857192993, 0.0169016718864441, 0.0178416967391968, 0.0124849081039429, 0.0122348666191101, 0.00828582048416138,0]]) mix.kij_cubic(Kij) eos = preos(mix, 'qmr')

n=2 T=123.396045640827 x1=np.linspace(0,1,n) x2=np.linspace(0,1,n) x3=np.linspace(0,1,n) x4=np.linspace(0,1,n) x5=np.linspace(0,1,n) x6=np.linspace(0,1,n) x7=np.linspace(0,1,n) x8=np.linspace(0,1,n) X=np.array([x1,x2,x3,x4,x5,x6,x7,x8,1-x1+x2+x3+x4+x5+x6+x7+x8]) Y=np.zeros([9,n]) P=np.zeros([n]) i=0 P0=1 #chute inicial de Temperatura Y0=np.array([0.883983915077912,0.053542087690542,0.0325279486946615,0.0135848429498525,0.00786401076896055,0.00273475390472645,0.00244812590645285,0.00154630024103214,0.00176801476585997]) #Chute inicial de composição Y[:,i],P[i] = bubblePy(Y0,P0,X[:,i],T,eos)

for i in range (1, n): Y[:,i], P[i] = bubblePy(Y[:,i-1], P[i-1], X[:, i], T, eos)
print (x1,P)

[Danlira1991]

I've been trying without sucess is write a code that generate a full PxT phase diagram using BubblePy and DewPx for this system (c1,c2,c3,n-c4,ic4,n-c5,i-c5,c6,c7+) in oder to compare with Aspen Hysys solution. I have little knowledge of the Python language so making the union between bubblePy and dewPx and generating the complete diagram has been challenging for me.

[Danlira1991]

Hi Gustavo, I redid the case I shared and fixed some input errors. The simulations were run using bubblePy for the nine component system. For the initial pressure kicks, I used those calculated by hysys, as if they were experimental (or reference) data. The first initial guess for the composition of the vapor phase was also taken from Hysys and for the others, the compositions calculated previously were used successively. The results obtained are shown in the figure below. In the cricondentherm region the results start to diverge from the Hysys solution and from this point on, no bubble point can be calculated, showing the error 'Array must not contain infs or NaNs'. Do you have any suggestions to make so that the entire bubble point curve is calculated like in Hysys? What about an automated way to do these bubble point calculations like you did in the article 'Phasepy: A Python based framework...' in listing 12 for the nine component case?

Ps: Corrected set-up for the first point:

import numpy as np from phasepy import mixture, component, preos, cubic from phasepy.equilibrium import bubblePy, dewPx, bubbleTy from phasepy.equilibrium import tpd_min, tpd_minimas

methane = component(name='methane', Tc=190.7, Pc=46.41, Zc=0.286, w=0.0115) ethane = component(name='ethane', Tc=305.43, Pc=48.84, Zc=0.279, w=0.0986) propane = component(name='propane', Tc=369.9, Pc=42.57, Zc=0.276, w=0.1524) i_butane = component(name='i_butane', Tc=408.05, Pc=36.48, Zc=0.282, w=0.1848) n_butane = component(name='n_butane', Tc=425.15, Pc=37.97, Zc=0.274, w=0.2010) i_pentane = component(name='i_pentane', Tc=460.35, Pc=33.34, w=0.2222) n_pentane = component(name='n_pentane', Tc=469.65, Pc=33.75, Zc=0.27, w=0.2539) n_Hexane = component(name='n_Hexane', Tc=507.85, Pc=30.32, Zc=0.266, w=0.3007) n7_plus = component(name='n7_plus', Tc=712.85, Pc=15.89, w=0.6530)

mix = methane + ethane+ propane +i_butane + n_butane+ i_pentane+n_pentane+n_Hexane+n7_plus Kij=np.array([[0, 0.00224137306213379, 0.00682878494262695, 0.0131134390830994, 0.0123046040534973, 0.017627477645874, 0.0179253816604614 ,0.0234740376472473, 0.0576416254043579], [0.00224137306213379, 0, 0.00125795602798462, 0.00457346439361572, 0.0040963888168335, 0.00741332769393921, 0.00760942697525024, 0.0114138126373291, 0.0382192730903625], [0.00682878494262695, 0.00125795602798462, 0, 0.00104051828384399, 0.000818967819213867, 0.00258338451385498, 0.00270050764083862, 0.0051419734954834 ,0.0260568857192993], [0.0131134390830994, 0.00457346439361572, 0.00104051828384399, 0, 0.0000133514404296875 , 0.000346183776855469, 0.000390052795410156, 0.00156527757644653, 0.0169016718864441], [0.0123046040534973, 0.0040963888168335, 0.000818967819213867, 0.0000133514404296875, 0, 0.00049513578414917, 0.000547230243682861, 0.00186634063720703, 0.0178416967391968],
[0.017627477645874, 0.00741332769393921, 0.00258338451385498, 0.000346183776855469, 0.00049513578414917,0, 0.0000012516975402832, 0.000439941883087158, 0.0124849081039429], [0.0179253816604614, 0.00760942697525024, 0.00270050764083862, 0.000390052795410156, 0.000547230243682861, 0.0000012516975402832, 0, 0.000393450260162354, 0.0122348666191101], [0.0234740376472473, 0.0114138126373291, 0.0051419734954834, 0.00156527757644653, 0.00186634063720703, 0.000439941883087158, 0.000393450260162354, 0,0.00828582048416138], [0.0576416254043579, 0.0382192730903625, 0.0260568857192993, 0.0169016718864441, 0.0178416967391968, 0.0124849081039429, 0.0122348666191101, 0.00828582048416138,0]]) mix.kij_cubic(Kij) eos = preos(mix, 'qmr')

x = np.array([0.42, 0.05, 0.05, 0.03, 0.02,0.01, 0.01, 0.1, 0.4]) T =123.396045640827 y0 = np.array([0,883983915562479, 5,35420872543362e-002, 3,25279484826942e-002, 1,35848429136660e-002, 7,86401075870640e-003, 2,73475391674585e-003, 2,44812592072292e-003, 1,54630026265731e-003, 1,76801492799193e-03]) P0 =1.888875172 bubblePy(y0, P0, x, T, eos, full_output=True)

[gustavochm]

Hey Danlira,

Sorry about my late response, I've been busy 🥲. I have to say that I am not an expert on generating saturation envelopes. For what I see the functions in phasepy are working at the initial conditions but the fail afterwards. I think this might be related on how the phase envelope is computed, in the book from Michelsen and Mollerrup (Chapter 12) they gave further details on how to construct the phase envelope.

They explain that when constructing the phase envelope sometimes they change for the calculation at a given temperature to pressure depending and viceversa depending on some sensitives equations. I don't have an automated implementation for that yet 😢

Btw, I'm actually taking the course where they teach these stuff but I don't know when I will have the time to actually implement it into phasepy.

Gustavo

[Danlira1991]

Hey gustavochm, Thank you for the recomendation! It really helps.

I try to use good_initial and all points could be calculated. But I could not full understand what does in code when I choose good_initial. Can you help me on this?

Thanks once again!

[gustavochm]

Hey Danlira,

So basically, the bubble and dew points functions included in phasepy include two methods to solve the phase equilibria.

1. The first method uses an Newton method to update the temperature or pressure in a outer loop while the compositions are internally updated are using successive substitution. This method is of slow convergence rate, it is reasonable good at low pressure and when -poor- initial guesses are available.
2. The second method solves a mutidimensional systems of equations (see image bellow). This system of equation was proposed by Michelsen to compute the phase envelope. It is usually of a higher convergence rate but it requires a good initial guess to converge.

In phase the strategy is the following. If good_initial=False the function will first try with method 1. If no solution has been found after a given number of iterations the function will switch to method 2. If good_initial=True then the function will automatically go for method 2.

I hope that makes it easier to understand was it is going on under the hood.

Gustavo

[Danlira1991]

Hey Danlira,

So basically, the bubble and dew points functions included in phasepy include two methods to solve the phase equilibria.

1. The first method uses an Newton method to update the temperature or pressure in a outer loop while the compositions are internally updated are using successive substitution. This method is of slow convergence rate, it is reasonable good at low pressure and when -poor- initial guesses are available.
2. The second method solves a mutidimensional systems of equations (see image bellow). This system of equation was proposed by Michelsen to compute the phase envelope. It is usually of a higher convergence rate but it requires a good initial guess to converge.

In phase the strategy is the following. If good_initial=False the function will first try with method 1. If no solution has been found after a given number of iterations the function will switch to method 2. If good_initial=True then the function will automatically go for method 2.

I hope that makes it easier to understand was it is going on under the hood.

Gustavo

Yeah! Thank you, thats what I was thinking.

[jhalzate]

I'm working on my own phase envelope printer. Will make the code available later after I test it properly and package it...

[gustavochm]

I'm working on my own phase envelope printer. Will make the code available later after I test it properly and package it...

That looks super good! Glad you are going to share it :)!