using Pkg
Pkg.activate("..")
using Clapeyron, PyCall, CSV, Tables
# using Pkg
# Pkg.add("PyPlot")
import PyPlot; const plt = PyPlot
Activating project at `~/Library/CloudStorage/OneDrive-ImperialCollegeLondon/University/UROP/SAFT_codes/Clapeyron`
PyPlot
In this notebook, we will be giving examples of the various equilibrium and critical properties one can obtain in Clapeyron.jl
. As most of the previous notebooks have shown the more-classic $p-xy$ and $T-xy$ diagrams, this notebook will focus on the uncommon, more-complex types of diagrams.
As a first example, we consider a $pT$ isopleth for a four-component systems of a liquid-natural gas system.
Generating the models and defining the mole fraction for the isopleth:
model1 = GERG2008(["methane","butane","isobutane","pentane"])
model2 = PR(["methane","butane","isobutane","pentane"])
model3 = PCSAFT(["methane","butane","isobutane","pentane"])
model4 = SAFTgammaMie(["methane","butane","isobutane","pentane"])
models = [model1,model2,model3,model4]
x = [0.96611,0.01475,0.01527,0.00385];
Obtaining the dew and bubble points of the isopleth. Here, we guestimate the critical point of the mixture and re-use the initial guess to reduce the computational time:
p = []
T = []
Tcrit = [188.039,204.173,202.694,207.7]
Tmax = [267,265.5,267,267]
for j=1:4
T1 = range(150,Tcrit[j],length=400)
T2 = range(200,Tmax[j],length=400)
T3 = range(Tmax[j],215,length=400)
bub1 = []
bub2 = []
bub3 = []
v1 = []
v2 = []
for i ∈ 1:400
if i ==1
append!(bub1,[bubble_pressure(models[j],T1[i],x)])
v1 = [log10(bub1[i][2]),log10(bub1[i][3]),bub1[i][4][1],bub1[i][4][2],bub1[i][4][3],bub1[i][4][4]]
append!(bub2,[dew_pressure(models[j],T2[i],x)])
v2 = [log10(bub2[i][2]),log10(bub2[i][3]),bub2[i][4][1],bub2[i][4][2],bub2[i][4][3],bub2[i][4][4]]
else
append!(bub1,[bubble_pressure(models[j],T1[i],x;v0=v1)])
v1 = [log10(bub1[i][2]),log10(bub1[i][3]),bub1[i][4][1],bub1[i][4][2],bub1[i][4][3],bub1[i][4][4]]
append!(bub2,[dew_pressure(models[j],T2[i],x;v0=v2)])
v2 = [log10(bub2[i][2]),log10(bub2[i][3]),bub2[i][4][1],bub2[i][4][2],bub2[i][4][3],bub2[i][4][4]]
end
end
v3 = deepcopy(v2)
v3[2] = 1.1*v3[2]
for i ∈ 1:400
#
if i ==1
append!(bub3,[dew_pressure(models[j],T3[i],x;v0=v3)])
v3 = [log10(bub3[i][2]),log10(bub3[i][3]),bub3[i][4][1],bub3[i][4][2],bub3[i][4][3],bub3[i][4][4]]
else
append!(bub3,[dew_pressure(models[j],T3[i],x;v0=v3)])
v3 = [log10(bub3[i][2]),log10(bub3[i][3]),bub3[i][4][1],bub3[i][4][2],bub3[i][4][3],bub3[i][4][4]]
end
end
p1 = [bub1[i][1] for i ∈ 1:400]
p2 = [bub2[i][1] for i ∈ 1:400]
p3 = [bub3[i][1] for i ∈ 1:400]
A = []
append!(A,p1)
append!(A,reverse(p3))
append!(A,reverse(p2))
append!(p,[A])
B = []
append!(B,T1)
append!(B,reverse(T3))
append!(B,reverse(T2))
append!(T,[B])
end
Note that we've split the isopleth into three sections: 2 dew curves and 1 bubble curve. This is primarily because, at a given temperature and composition, there are two possible dew points for this system. As such, we must scan the 'lower' dew curve first before switching the direction of the scan to obtain the 'upper' dew curve.
Plotting:
plt.clf()
plt.plot(T[1],p[1]./1e6,label="GERG-2008",linestyle="--")
plt.plot(T[2],p[2]./1e6,label="PR",linestyle=":")
plt.plot(T[3],p[3]./1e6,label="PC-SAFT",linestyle="-.")
plt.plot(T[4],p[4]./1e6,label="SAFT-γ Mie",linestyle=(0, (3, 1, 1, 1)))
plt.legend(loc="upper left",frameon=false,fontsize=12)
plt.xlabel("Temperature / K",fontsize=16)
plt.ylabel("Pressure / MPa",fontsize=16)
plt.ylim([0.,10.])
plt.xlim([150,300])
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
display(plt.gcf())
Unsurprisingly, GERG-2008 performs the best in this system; even more-surprising is how similar the predictions from the other equations of state are.
Ramping up the difficulty, we now examine a 12 component liquid natural gas system. This is to highlight that Clapeyron.jl
is able to handle systems of particularly large size. For some of these equations of state, these are the largest systems they've ever been shown to model.
Generating the models and setting the mole fraction for the isopleth:
model1 = GERG2008(["methane","nitrogen","carbon dioxide","ethane","propane","butane","isobutane",
"pentane","isopentane","hexane","heptane","octane"])
model2 = SRK(["methane","nitrogen","carbon dioxide","ethane","propane","butane","isobutane",
"pentane","isopentane","hexane","heptane","octane"])
model3 = PCSAFT(["methane","nitrogen","carbon dioxide","ethane","propane","butane","isobutane",
"pentane","isopentane","hexane","heptane","octane"])
model4 = SAFTgammaMie(["methane","nitrogen","carbon dioxide","ethane","propane","butane","isobutane",
"pentane","isopentane","hexane","heptane","octane"])
models = [model1,model2,model3,model4]
x = [0.833482,0.056510,0.002840,0.075260,0.020090,0.005200,0.003050,0.001440,0.001200,0.000680,0.000138,0.000110];
As before, we scan the isopleth in three chunks; again, we guestimate the critical point and re-use initial guesses:
p = []
T = []
Tcrit = [214,214,214,218]
Tmax = [272.5,272.5,273,273]
for j=1:4
# crit = crit_mix(models[j],x)
T1 = range(150,Tcrit[j],length=100)
T2 = range(225,Tmax[j],length=100)
T3 = range(Tmax[j],220,length=100)
bub1 = []
bub2 = []
bub3 = []
v1 = []
v2 = []
for i ∈ 1:100
if i ==1
append!(bub1,[bubble_pressure(models[j],T1[i],x)])
v1 = [log10(bub1[i][2]),log10(bub1[i][3])]
append!(v1,bub1[i][4])
append!(bub2,[dew_pressure(models[j],T2[i],x)])
v2 = [log10(bub2[i][2]),log10(bub2[i][3])]
append!(v2,bub2[i][4])
else
append!(bub1,[bubble_pressure(models[j],T1[i],x;v0=v1)])
v1 = [log10(bub1[i][2]),log10(bub1[i][3])]
append!(v1,bub1[i][4])
append!(bub2,[dew_pressure(models[j],T2[i],x;v0=v2)])
v2 = [log10(bub2[i][2]),log10(bub2[i][3])]
append!(v2,bub2[i][4])
end
end
v3 = deepcopy(v2)
v3[2] = 1.1*v3[2]
for i ∈ 1:100
if i ==1
append!(bub3,[dew_pressure(models[j],T3[i],x;v0=v3)])
v3 = [log10(bub3[i][2]),log10(bub3[i][3])]
append!(v3,bub3[i][4])
else
append!(bub3,[dew_pressure(models[j],T3[i],x;v0=v3)])
v3 = [log10(bub3[i][2]),log10(bub3[i][3])]
append!(v3,bub3[i][4])
end
end
p1 = [bub1[i][1] for i ∈ 1:100]
p2 = [bub2[i][1] for i ∈ 1:100]
p3 = [bub3[i][1] for i ∈ 1:100]
A = []
append!(A,p1)
append!(A,reverse(p3))
append!(A,reverse(p2))
append!(p,[A])
B = []
append!(B,T1)
append!(B,reverse(T3))
append!(B,reverse(T2))
append!(T,[B])
end
Plotting:
plt.clf()
plt.plot(T[1],p[1]./1e6,label="GERG-2008",linestyle="--")
plt.plot(T[2],p[2]./1e6,label="SRK",linestyle=":")
plt.plot(vcat(T[3][1:99],T[3][102:end]),vcat(p[3][1:99],p[3][102:end])./1e6,label="PC-SAFT",linestyle="-.")
plt.plot(T[4],p[4]./1e6,label="SAFT-γ Mie",linestyle=(0, (3, 1, 1, 1)))
plt.legend(loc="upper left",frameon=false,fontsize=12)
plt.xlabel("Temperature / K",fontsize=16)
plt.ylabel("Pressure / MPa",fontsize=16)
plt.ylim([0.,10.])
plt.xlim([150,300])
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
display(plt.gcf())
We've shown in previous notebooks how Clapeyron.jl
can plot the entire phase diagram for a single component; in a binary component the `complete' phase diagram is usually referred to as the $pT$ projection where we include the saturation curve for each component and the critical point of the mixture between these two curves. We first choose to model the carbon dioxide + ethane system as the critical curve here has an unusual shape.
Generating the models:
model1 = SRK(["carbon dioxide","ethane"])
model2 = PSRK(["carbon dioxide","ethane"])
model3 = GERG2008(["carbon dioxide","ethane"])
model4 = SAFTVRMie(["carbon dioxide","ethane"])
model5 = sPCSAFT(["carbon dioxide","ethane"])
models = [model1,model2,model3,model4,model5];
We first obtain the saturation curve for each component using each equation of state:
p_1 = []
p_2 = []
T_1 = []
T_2 = []
for j=1:5
pure = split_model(models[j])
crit = crit_pure.(pure)
T = []
sat = []
for i=1:2
append!(T,[range(0.7*crit[i][1],crit[i][1],length=100)])
append!(sat,[saturation_pressure.(pure[i],T[i])])
end
append!(p_1,[[sat[1][i][1] for i ∈ 1:100]])
append!(p_2,[[sat[2][i][1] for i ∈ 1:100]])
append!(T_1,[T[1]])
append!(T_2,[T[2]])
end
We next solve for the critical points of the mixture using crit_mix
where one specifies the mole fraction of the mixture:
x = range(1e-5,1-1e-5,length=100)
X = Clapeyron.FractionVector.(x)
T_crit = []
p_crit = []
for j=1:5
T = []
p = []
v0 = []
for i=1:100
if i==1
mix_crit = crit_mix(models[j],X[i])
v0 = [log10(mix_crit[3]),mix_crit[1]]
append!(T,mix_crit[1])
append!(p,mix_crit[2])
else
mix_crit = crit_mix(models[j],X[i];v0=v0)
v0 = [log10(mix_crit[3]),mix_crit[1]]
append!(T,mix_crit[1])
append!(p,mix_crit[2])
end
end
append!(T_crit,[T])
append!(p_crit,[p])
end
To reduce computational time, we re-use initial guesses.
Plotting:
plt.clf()
plt.plot(T_1[1],p_1[1]./1e6,label="",linestyle="-",color="r")
plt.plot(T_1[2],p_1[2]./1e6,label="",linestyle="-",color="b")
plt.plot(T_1[3],p_1[3]./1e6,label="",linestyle="-",color="g")
plt.plot(T_1[4],p_1[4]./1e6,label="",linestyle="-",color="purple")
plt.plot(T_1[5],p_1[5]./1e6,label="",linestyle="-",color="darkblue")
plt.plot(T_2[1],p_2[1]./1e6,label="",linestyle="-",color="r")
plt.plot(T_2[2],p_2[2]./1e6,label="",linestyle="-",color="b")
plt.plot(T_2[3],p_2[3]./1e6,label="",linestyle="-",color="g")
plt.plot(T_2[4],p_2[4]./1e6,label="",linestyle="-",color="purple")
plt.plot(T_2[5],p_2[5]./1e6,label="",linestyle="-",color="darkblue")
plt.plot(T_1[1][end],p_1[1][end]./1e6,label="",marker="^",linestyle="",color="r")
plt.plot(T_1[2][end],p_1[2][end]./1e6,label="",marker="^",linestyle="",color="b")
plt.plot(T_1[3][end],p_1[3][end]./1e6,label="",marker="^",linestyle="",color="g")
plt.plot(T_1[4][end],p_1[4][end]./1e6,label="",marker="^",linestyle="",color="purple")
plt.plot(T_1[5][end],p_1[5][end]./1e6,label="",marker="^",linestyle="",color="darkblue")
plt.plot(T_2[1][end],p_2[1][end]./1e6,label="",marker="^",linestyle="",color="r")
plt.plot(T_2[2][end],p_2[2][end]./1e6,label="",marker="^",linestyle="",color="b")
plt.plot(T_2[3][end],p_2[3][end]./1e6,label="",marker="^",linestyle="",color="g")
plt.plot(T_2[4][end],p_2[4][end]./1e6,label="",marker="^",linestyle="",color="purple")
plt.plot(T_2[5][end],p_2[5][end]./1e6,label="",marker="^",linestyle="",color="darkblue")
plt.plot(T_crit[1],p_crit[1]./1e6,label="SRK",linestyle=":",color="r")
plt.plot(T_crit[2],p_crit[2]./1e6,label="PSRK",linestyle="--",color="b")
plt.plot(T_crit[3],p_crit[3]./1e6,label="GERG-2008",linestyle="-.",color="g")
plt.plot(T_crit[5],p_crit[5]./1e6,label="sPC-SAFT",linestyle=(0, (5, 1)),color="darkblue")
plt.plot(T_crit[4],p_crit[4]./1e6,label="SAFT-VR Mie",linestyle=(0, (3, 1, 1, 1)),color="purple")
plt.legend(loc="upper left",frameon=false,fontsize=12)
plt.xlabel("Temperature / K",fontsize=16)
plt.ylabel("Pressure / MPa",fontsize=16)
plt.ylim([3.,9.])
plt.xlim([270,320])
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
display(plt.gcf())
In the previous case, we had two components with very similar critical points; now we examine the inverse case where the critical points are different.
Generating the models:
model1 = VTPR(["ethane","decane"])
model2 = GERG2008(["ethane","decane"])
model3 = SAFTVRMie(["ethane","decane"])
model4 = SAFTgammaMie(["ethane","decane"])
model5 = softSAFT(["ethane","decane"])
models = [model1,model2,model3,model4,model5];
Obtaining the saturation curve for the pure components:
p_1 = []
p_2 = []
T_1 = []
T_2 = []
for j=1:5
pure = split_model(models[j])
crit = crit_pure.(pure)
T = []
sat = []
for i=1:2
append!(T,[range(0.7*crit[i][1],crit[i][1],length=100)])
append!(sat,[saturation_pressure.(pure[i],T[i])])
end
append!(p_1,[[sat[1][i][1] for i ∈ 1:100]])
append!(p_2,[[sat[2][i][1] for i ∈ 1:100]])
append!(T_1,[T[1]])
append!(T_2,[T[2]])
end
Obtaining the critical point of the mixture at different compositions, re-using initial guesses:
x = 1 .-range(1e-5,1-1e-4,length=200)
X = Clapeyron.FractionVector.(x)
T_crit = []
p_crit = []
for j=1:5
T = []
p = []
v0 = []
for i=1:200
if i==1
mix_crit = crit_mix(models[j],X[i])
v0 = [log10(mix_crit[3]),mix_crit[1]]
append!(T,mix_crit[1])
append!(p,mix_crit[2])
else
mix_crit = crit_mix(models[j],X[i];v0=v0)
v0 = [log10(mix_crit[3]),mix_crit[1]]
append!(T,mix_crit[1])
append!(p,mix_crit[2])
end
end
append!(T_crit,[T])
append!(p_crit,[p])
end
Plotting:
plt.clf()
plt.plot(T_1[1],p_1[1]./1e6,label="",linestyle="-",color="r")
plt.plot(T_1[2],p_1[2]./1e6,label="",linestyle="-",color="b")
plt.plot(T_1[3],p_1[3]./1e6,label="",linestyle="-",color="g")
plt.plot(T_1[4],p_1[4]./1e6,label="",linestyle="-",color="purple")
plt.plot(T_1[5],p_1[5]./1e6,label="",linestyle="-",color="lightblue")
plt.plot(T_2[1],p_2[1]./1e6,label="",linestyle="-",color="r")
plt.plot(T_2[2],p_2[2]./1e6,label="",linestyle="-",color="b")
plt.plot(T_2[3],p_2[3]./1e6,label="",linestyle="-",color="g")
plt.plot(T_2[4],p_2[4]./1e6,label="",linestyle="-",color="purple")
plt.plot(T_2[5],p_2[5]./1e6,label="",linestyle="-",color="lightblue")
plt.plot(T_1[1][end],p_1[1][end]./1e6,label="",marker="^",linestyle="",color="r")
plt.plot(T_1[2][end],p_1[2][end]./1e6,label="",marker="^",linestyle="",color="b")
plt.plot(T_1[3][end],p_1[3][end]./1e6,label="",marker="^",linestyle="",color="g")
plt.plot(T_1[4][end],p_1[4][end]./1e6,label="",marker="^",linestyle="",color="purple")
plt.plot(T_1[5][end],p_1[5][end]./1e6,label="",marker="^",linestyle="",color="lightblue")
plt.plot(T_2[1][end],p_2[1][end]./1e6,label="",marker="^",linestyle="",color="r")
plt.plot(T_2[2][end],p_2[2][end]./1e6,label="",marker="^",linestyle="",color="b")
plt.plot(T_2[3][end],p_2[3][end]./1e6,label="",marker="^",linestyle="",color="g")
plt.plot(T_2[4][end],p_2[4][end]./1e6,label="",marker="^",linestyle="",color="purple")
plt.plot(T_2[5][end],p_2[5][end]./1e6,label="",marker="^",linestyle="",color="lightblue")
plt.plot(T_crit[1],p_crit[1]./1e6,label="VTPR",linestyle=":",color="r")
plt.plot(T_crit[2],p_crit[2]./1e6,label="GERG-2008",linestyle="--",color="b")
plt.plot(T_crit[5],p_crit[5]./1e6,label="soft-SAFT",linestyle=(0, (5, 1)),color="lightblue")
plt.plot(T_crit[3],p_crit[3]./1e6,label="SAFT-VR Mie",linestyle="-.",color="g")
plt.plot(T_crit[4],p_crit[4]./1e6,label="SAFT-γ Mie",linestyle=(0, (3, 1, 1, 1)),color="purple")
plt.legend(loc="upper right",frameon=false,fontsize=12)
plt.xlabel("Temperature / K",fontsize=16)
plt.ylabel("Pressure / MPa",fontsize=16)
plt.ylim([0.,20.])
plt.xlim([280,700])
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
display(plt.gcf())
In the previous examples, we only considered type-I mixtures. To illustrate how Clapeyron.jl
can indeed handle different binary phase types, we consider the methanol + cyclohexane system implemented in PC-SAFT. Here, we first obtain the upper critical end point of the mixture (where the UCST curve meets the VLLE and azeotrope curve). From here, we scan for the UCST, VLLE and azeotrope curves. The saturation and critical curves are obtained as before:
model = PCSAFT(["methanol","cyclohexane"])
pure = Clapeyron.split_model(model)
crit = crit_pure.(pure);
UCEP = UCEP_mix(model);
T_1 = range(310,crit[1][1],length=200)
T_2 = range(310,crit[2][1],length=200)
T_3 = range(310,UCEP[1]*0.999,length=200)
T_4 = []
T_5 = range(UCEP[1],520,length=200)
T_6 = range(UCEP[1],UCEP[1]*1.45,length=200)
x = range(1e-30,1-1e-6,length=200)
X = Clapeyron.FractionVector.(x)
p_1 = []
p_2 = []
p_3 = []
p_4 = []
p_5 = []
p_6 = []
v0_1 = 0.
v0_2 = 0.
v0_3 = 0.
v0_4 = 0.
v0_5 = 0.
v0_6 = 0.
v_l_5 = []
v_v_5 = []
x_5 = []
for i in 1:200
if i == 1
sat_1 = saturation_pressure(pure[1],T_1[i])
v0_1 = log10.([sat_1[2],sat_1[3]])
append!(p_1,sat_1[1])
sat_2 = saturation_pressure(pure[2],T_2[i])
v0_2 = log10.([sat_2[2],sat_2[3]])
append!(p_2,sat_2[1])
three = VLLE_pressure(model,T_3[i])
v0_3 = [log10(three[2]),log10(three[3]),log10(three[4]),three[5][1],three[6][1],three[7][1]]
append!(p_3,three[1])
crit_1 = crit_mix(model,X[i])
v0_4 = [log10(crit_1[3]),crit_1[1]]
append!(p_4,crit_1[2])
append!(T_4,crit_1[1])
az = azeotrope_pressure(model,T_5[i])
v0_5 = az[4]
append!(p_5,az[1])
append!(v_l_5,az[2])
append!(v_v_5,az[3])
append!(x_5,az[4][1])
UCST_1 = UCST_mix(model,T_6[i];v0=[log10(UCEP[3]),UCEP[5]])
v0_6 = [log10(UCST_1[2]),UCST_1[3]]
append!(p_6,UCST_1[1])
else
sat_1 = saturation_pressure(pure[1],T_1[i],v0_1)
v0_1 = log10.([sat_1[2],sat_1[3]])
append!(p_1,sat_1[1])
sat_2 = saturation_pressure(pure[2],T_2[i],v0_2)
v0_2 = log10.([sat_2[2],sat_2[3]])
append!(p_2,sat_2[1])
three = VLLE_pressure(model,T_3[i];v0=v0_3)
v0_3 = [log10(three[2]),log10(three[3]),log10(three[4]),three[5][1],three[6][1],three[7][1]]
append!(p_3,three[1])
crit_1 = crit_mix(model,X[i];v0=v0_4)
v0_4 = [log10(crit_1[3]),crit_1[1]]
append!(p_4,crit_1[2])
append!(T_4,crit_1[1])
az = azeotrope_pressure(model,T_5[i];v0=v0_5)
v0_5 = az[4]
append!(p_5,az[1])
append!(v_l_5,az[2])
append!(v_v_5,az[3])
append!(x_5,az[4][1])
UCST_1 = UCST_mix(model,T_6[i];v0=v0_6)
v0_6 = [log10(UCST_1[2]),UCST_1[3]]
append!(p_6,UCST_1[1])
end
end
Plotting:
plt.clf()
plt.semilogy(crit[1][1],crit[1][2]./1e6,label="Critical (pure)",linestyle="",color="red",marker="^")
plt.plot(crit[2][1],crit[2][2]./1e6,label="",linestyle="",color="red",marker="^")
plt.plot(UCEP[1],UCEP[2]./1e6,label="UCEP",linestyle="",color="k",marker="s")
plt.plot(T_1,p_1./1e6,label="Saturation",color="k")
plt.plot(T_2,p_2./1e6,label="",color="k")
plt.plot(T_4,p_4./1e6,label="Critical (mixture)",color="r")
plt.plot(T_3,p_3./1e6,label="VLLE",color="blue")
plt.plot(T_6,p_6./1e6,label="UCST",color="green")
plt.plot(T_5,p_5./1e6,label="Azeotrope",color="purple")
plt.legend(loc="lower right",frameon=false,fontsize=11)
plt.xlabel("Temperature / K",fontsize=16)
plt.ylabel("Pressure / MPa",fontsize=16)
plt.ylim([3e-3,1e3])
plt.xlim([310,600])
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
display(plt.gcf())
In all the previous notebooks, we only considered VLE envelopes of binary components. However, binary phase diagrams can become more-complicated than this as we can include VLE, LLE and VLLE. To illustrate that these diagrams can be plotted in Clapeyron.jl
, we consider the methanol + ethane system.
Generating the models:
model1 = PCSAFT(["methanol","ethane"])
model2 = SAFTgammaMie(["methanol","ethane"])
model3 = PR(["methanol","ethane"])
PR{BasicIdeal, PRAlpha, NoTranslation, vdW1fRule} with 2 components: "methanol" "ethane" Contains parameters: a, b, Tc, Pc, Mw
By providing a good initial guess, we can obtain the VLLE point, from which we can scan for the VLE and LLE envelopes (as before, we re-use initial guesses). For PC-SAFT:
VLLE_PC=VLLE_pressure(model1,298.15;v0 = [log10(4.8654636008937104e-5),log10(7.709740609838804e-5),log10(0.00035501867564056083),0.64,0.08,0.007])
x1_PC = range(1e-6,VLLE_PC[6][1],length=100)
x2_PC = range(VLLE_PC[5][1],1-1e-6,length=100)
x3_PC = range(VLLE_PC[6][1],5*VLLE_PC[6][1],length=100)
X1 = Clapeyron.FractionVector.(x1_PC)
X2 = Clapeyron.FractionVector.(x2_PC)
X3 = Clapeyron.FractionVector.(x3_PC)
bub1 = bubble_pressure.(model1,298.15,X1)
bub2 = bubble_pressure.(model1,298.15,X2)
LLE = LLE_pressure.(model1,298.15,X3)
p1_PC = [bub1[i][1] for i ∈ 1:100]
p2_PC = [bub2[i][1] for i ∈ 1:100]
p3_PC = [LLE[i][1] for i ∈ 1:100]
y1_PC = [bub1[i][4][1] for i ∈ 1:100]
y2_PC = [bub2[i][4][1] for i ∈ 1:100]
y3_PC = [LLE[i][4][1] for i ∈ 1:100];
For SAFT-γ Mie:
VLLE_gamma=VLLE_pressure(model2,298.15;v0 = [log10(4.8654636008937104e-5),log10(8.709740609838804e-5),log10(0.00035501867564056083),0.64,0.08,0.007])
x1_gamma = range(1e-6,VLLE_gamma[6][1],length=100)
x2_gamma = range(VLLE_gamma[5][1],1-1e-6,length=100)
x3_gamma = range(VLLE_gamma[6][1],1.2*VLLE_gamma[6][1],length=100)
X1 = Clapeyron.FractionVector.(x1_gamma)
X2 = Clapeyron.FractionVector.(x2_gamma)
X3 = Clapeyron.FractionVector.(x3_gamma)
bub1 = bubble_pressure.(model2,298.15,X1)
bub2 = bubble_pressure.(model2,298.15,X2)
LLE = LLE_pressure.(model2,298.15,X3)
p1_gamma = [bub1[i][1] for i ∈ 1:100]
p2_gamma = [bub2[i][1] for i ∈ 1:100]
p3_gamma = [LLE[i][1] for i ∈ 1:100]
y1_gamma = [bub1[i][4][1] for i ∈ 1:100]
y2_gamma = [bub2[i][4][1] for i ∈ 1:100]
y3_gamma = [LLE[i][4][1] for i ∈ 1:100];
For PR:
VLLE_PR=VLLE_pressure(model3,298.15;v0 = [log10(5.0829874375196166e-5),log10(6.433971553612239e-5),log10(0.00043337290431935867),0.64,0.20,0.01])
x1_PR = range(1e-6,VLLE_PR[6][1],length=100)
x2_PR = range(VLLE_PR[5][1],1-1e-6,length=100)
x3_PR = range(VLLE_PR[6][1],2*VLLE_PR[6][1],length=100)
X1 = Clapeyron.FractionVector.(x1_PR)
X2 = Clapeyron.FractionVector.(x2_PR)
X3 = Clapeyron.FractionVector.(x3_PR)
v1 = [log10(VLLE_PR[3]),log10(VLLE_PR[4]),VLLE_PR[7][1],VLLE_PR[7][2]]
v2 = [log10(VLLE_PR[2]),log10(VLLE_PR[4]),VLLE_PR[7][1],VLLE_PR[7][2]]
v3 = [log10(VLLE_PR[3]),log10(VLLE_PR[2]),VLLE_PR[5][1]]
p1_PR = []
p2_PR = []
p3_PR = []
y1_PR = []
y2_PR = []
y3_PR = []
for i=1:100
bub1 = bubble_pressure(model3,298.15,X1[i];v0=v1)
bub2 = bubble_pressure(model3,298.15,X2[i];v0=v2)
LLE = LLE_pressure(model3,298.15,X3[i];v0=v3)
v1 = [log10(bub1[2]),log10(bub1[3]),bub1[4][1],bub1[4][2]]
v2 = [log10(bub2[2]),log10(bub2[3]),bub2[4][1],bub2[4][2]]
v3 = [log10(LLE[2]),log10(LLE[3]),LLE[4][1]]
append!(p1_PR,bub1[1])
append!(p2_PR,bub2[1])
append!(p3_PR,LLE[1])
append!(y1_PR,bub1[4][1])
append!(y2_PR,bub2[4][1])
append!(y3_PR,LLE[4][1])
end
Plotting:
plt.clf()
plt.plot(y1_PR,p1_PR./1e6,label="PR",linestyle="-",color="g")
plt.plot(x1_PR,p1_PR./1e6,label="",linestyle="-",color="g")
plt.plot(y2_PR,p2_PR./1e6,label="",linestyle="-",color="g")
plt.plot(x2_PR,p2_PR./1e6,label="",linestyle="-",color="g")
plt.plot(y3_PR,p3_PR./1e6,label="",linestyle="-",color="g")
plt.plot(x3_PR,p3_PR./1e6,label="",linestyle="-",color="g")
plt.plot([VLLE_PR[5][1],VLLE_PR[6][1],VLLE_PR[7][1]],[VLLE_PR[1],VLLE_PR[1],VLLE_PR[1]]./1e6,label="",linestyle="-.",color="g")
plt.plot(y1_PC,p1_PC./1e6,label="PC-SAFT",linestyle="-",color="b")
plt.plot(x1_PC,p1_PC./1e6,label="",linestyle="-",color="b")
plt.plot(y2_PC,p2_PC./1e6,label="",linestyle="-",color="b")
plt.plot(x2_PC,p2_PC./1e6,label="",linestyle="-",color="b")
plt.plot(y3_PC,p3_PC./1e6,label="",linestyle="-",color="b")
plt.plot(x3_PC,p3_PC./1e6,label="",linestyle="-",color="b")
plt.plot([VLLE_PC[5][1],VLLE_PC[6][1],VLLE_PC[7][1]],[VLLE_PC[1],VLLE_PC[1],VLLE_PC[1]]./1e6,label="",linestyle="-.",color="b")
plt.plot(y1_gamma,p1_gamma./1e6,label="SAFT-γ Mie",linestyle="-",color="r")
plt.plot(x1_gamma,p1_gamma./1e6,label="",linestyle="-",color="r")
plt.plot(y2_gamma,p2_gamma./1e6,label="",linestyle="-",color="r")
plt.plot(x2_gamma,p2_gamma./1e6,label="",linestyle="-",color="r")
plt.plot(y3_gamma,p3_gamma./1e6,label="",linestyle="-",color="r")
plt.plot(x3_gamma,p3_gamma./1e6,label="",linestyle="-",color="r")
plt.plot([VLLE_gamma[5][1],VLLE_gamma[6][1],VLLE_gamma[7][1]],[VLLE_gamma[1],VLLE_gamma[1],VLLE_gamma[1]]./1e6,label="",linestyle="-.",color="r")
plt.legend(loc="upper right",frameon=false,fontsize=11)
plt.xlabel("composition of 1",fontsize=16)
plt.ylabel("Pressure / MPa",fontsize=16)
plt.ylim([0,6])
plt.xlim([0,1])
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
display(plt.gcf())
Showing how we can also produce $T-xxy$ diagrams which contain VLE, LLE and VLLE regions using Clapeyron.jl
, we consider the methanol + cyclohexane system.
Generating the models:
model1 = SAFTgammaMie(["methanol","cyclohexane"])
model2 = PCSAFT(["methanol","cyclohexane"])
model3 = SRK(["methanol","cyclohexane"])
RK{BasicIdeal, SoaveAlpha, NoTranslation, vdW1fRule} with 2 components: "methanol" "cyclohexane" Contains parameters: a, b, Tc, Pc, Mw
By providing a good initial guess, we can obtain the VLLE point, from which we can scan for the VLE and LLE envelopes (as before, we re-use initial guesses). For SAFT-γ Mie:
VLLE_gamma = VLLE_pressure(model1,325.8576;v0=[log10(4.366652011766123e-5),
log10(0.0001098063430027416),
log10(0.023701084645953704),
0.98,0.05,0.62])
x1_gamma = range(1e-5,VLLE_gamma[6][1],length=100)
x2_gamma = range(VLLE_gamma[5][1],1-1e-5,length=100)
x3_gamma = range(VLLE_gamma[5][1],VLLE_gamma[5][1]*1.008,length=100)
X1 = Clapeyron.FractionVector.(x1_gamma)
X2 = Clapeyron.FractionVector.(x2_gamma)
X3 = Clapeyron.FractionVector.(x3_gamma)
bub1 = []
bub2 = []
LLE = []
T0 = 325.8576
for i ∈ 1:100
append!(bub1, [bubble_temperature(model1,1.013e5,X1[i])])
append!(bub2, [bubble_temperature(model1,1.013e5,X2[i])])
append!(LLE ,[LLE_temperature(model1,1.013e5,X3[i])])
T0 = LLE[i][1]
end
T1_gamma = [bub1[i][1] for i ∈ 1:100]
T2_gamma = [bub2[i][1] for i ∈ 1:100]
T3_gamma = [LLE[i][1] for i ∈ 1:100]
y1_gamma = [bub1[i][4][1] for i ∈ 1:100]
y2_gamma = [bub2[i][4][1] for i ∈ 1:100]
y3_gamma = [LLE[i][4][1] for i ∈ 1:100];
For PC-SAFT:
x1_PC = range(1e-5,1-1e-5,length=200)
X = Clapeyron.FractionVector.(x)
x2 = range(0.05,0.4155,length=200)
X2 = Clapeyron.FractionVector.(x2)
T0 = 290
T1 = 340
LLE = []
bub = []
for i ∈ 1:200
append!(LLE,[LLE_temperature(model2,1.013e5,X2[i])])
T0 = LLE[i][1]
append!(bub, [bubble_temperature(model2,1.013e5,X[i])])
T1 = bub[i][1]
end
T1_PC = [bub[i][1] for i ∈ 1:200]
y1_PC = [bub[i][4][1] for i ∈ 1:200]
append!(y1_PC,reverse(x1_PC))
append!(T1_PC,reverse(T1_PC))
T2_PC = [LLE[i][1] for i ∈ 1:200]
x2_PC = [LLE[i][4][1] for i ∈ 1:200]
append!(x2_PC,reverse(x2))
append!(T2_PC,reverse(T2_PC));
For PR:
x1_PR = range(1e-5,1-1e-5,length=200)
X = Clapeyron.FractionVector.(x)
T0 = 352
bub = []
for i ∈ 1:200
append!(bub,[bubble_temperature(model3,1.013e5,X[i])])
end
T1_PR = [bub[i][1] for i ∈ 1:200]
y1_PR = [bub[i][4][1] for i ∈ 1:200]
append!(y1_PR,reverse(x1_PR))
append!(T1_PR,reverse(T1_PR));
Due to numerical errors, the LLE region was obtained using the tp_flash
method in Clapeyron.jl
:
T = range(290,317.3,length=100)
T = [T[i] for i ∈ 1:100]
method =
z = [0.65,0.35]
x1 = []
x2 = []
K0 = [4.5,0.125]
for i =1:100
flash = tp_flash(model3,1.013e5,T[i],z,RRTPFlash(K0=K0))
append!(x1,flash[1][1,1])
append!(x2,flash[1][2,1])
K0 = flash[1][2,:]./flash[1][1,:]
end
append!(x1,reverse(x2))
append!(T,reverse(T));
Plotting:
plt.clf()
plt.plot(y1_PR,T1_PR,label="SRK",color="g")
plt.plot(x1,T,label="",color="g")
plt.plot(y1_PC,T1_PC,label="PC-SAFT",color="b")
plt.plot(x2_PC,T2_PC,label="",color="b")
plt.plot(y1_gamma,T1_gamma,label="SAFT-γ Mie",color="r")
plt.plot(x1_gamma,T1_gamma,label="",color="r")
plt.plot(y2_gamma,T2_gamma,label="",color="r")
plt.plot(x2_gamma,T2_gamma,label="",color="r")
plt.plot(y3_gamma,T3_gamma,label="",color="r")
plt.plot(x3_gamma,T3_gamma,label="",color="r")
plt.plot(x3_gamma,T3_gamma,label="",color="r")
plt.plot([VLLE_gamma[5][1],VLLE_gamma[7][1],VLLE_gamma[6][1]],[325.8574,325.8574,325.8574],label="",linestyle="-.",color="r")
plt.legend(loc="upper right",frameon=false,fontsize=11)
plt.xlabel("composition of 1",fontsize=16)
plt.ylabel("Temperature / K",fontsize=16)
plt.xlim([0,1])
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.ylim(295,355)
display(plt.gcf())