Skip to content

Question: cubic_eos example doesn't work (?) #552

Description

@Emanukka

I was trying to reproduce the cubic eos example , but when I run it in my pc, it gives a different result in relation to what is in github notebook :

Image

(what should give):

Image

Also, I noted that in the same notebook example, the model2 is PatelTeja, which is different from PR:

Image

Otherwise, even though I change to PR with PR alpha function without volume translation, the result is the same.

Image

my code in vscode that i'm running:

using Pkg
# Pkg.activate("..")
using Clapeyron, PyCall
import PyPlot; const plt = PyPlot
# Generating the models
model1 = RK(["water"])
model2 = PR(["water"])
model3 = PR(["water"];alpha=TwuAlpha,translation=NoTranslation)
model4 = PR(["water"];alpha=TwuAlpha,translation=RackettTranslation)

# Concatenating them into one list
models = [model1,model2,model3,model4];

# Obtaining all of their critical points
crit = crit_pure.(models);

# Obtaining the saturation curve
N = 300
T = zeros(N,4)
p = zeros(N,4)
v_l = zeros(N,4)
v_v = zeros(N,4)
for i  1:4
    T[:,i] = range(285,crit[i][1],length=N)
    sat = saturation_pressure.(models[i],T[:,i])
    p[:,i] = [sat[i][1] for i  1:N]
    v_l[:,i] = [sat[i][2] for i  1:N]
    v_v[:,i] = [sat[i][3] for i  1:N]
end

T_exp = [300,320,340,360,380,400,420,440,460,480,500,520,540,560,580,600,620,640,647.096]

p_exp = [0.0035368,0.010546,0.027188,0.062194,0.12885,0.24577,0.4373,0.73367,1.1709,1.7905,2.6392,3.769,5.2369,7.1062,9.448,12.345,15.901,20.265,22.064]

ρ_l_exp = [55.315,54.919,54.371,53.698,52.918,52.038,51.064,49.994,48.824,47.545,46.145,44.603,42.889,40.956,38.725,36.048,32.577,26.729,17.873728]

ρ_v_exp = [0.0014204,0.0039778,0.0096808,0.021014,0.041537,0.076014,0.13055,0.21276,0.33209,0.50035,0.73265,1.0491,1.478,2.062,2.872,4.0434,5.9009,9.8331,17.873728];



# Plotting p-T curve
plt.clf()
plt.semilogy(T[:,1],p[:,1] ./1e6,label="RK",linestyle=":",color="r")
plt.semilogy(T[:,2],p[:,2] ./1e6,label="PR",linestyle="--",color="g")
plt.semilogy(T[:,3],p[:,3] ./1e6,label="PR{TwuAlpha}",linestyle="-.",color="b")
plt.semilogy(T[:,4],p[:,4] ./1e6,label="PR{TwuAlpha, RackettTranslation}",linestyle="-.")
plt.semilogy(T_exp,p_exp,label="Experimental",marker="o",linestyle="",color="k")
# Create a legend of all the existing plots using their labels as names
plt.legend(loc="lower right",frameon=false,fontsize=12) 
plt.xlabel("Temperature / K",fontsize=16)
plt.ylabel("Pressure / MPa",fontsize=16)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.xlim([300,650])
plt.ylim([5e-3,3e1])
display(plt.gcf())


# Plotting rho-T envelope
plt.clf()
plt.semilogx(1e-3 ./v_l[:,1],T[:,1],label="RK",linestyle=":",color="r")
plt.semilogx(1e-3 ./v_v[:,1],T[:,1],label="",linestyle=":",color="r")
plt.semilogx(1e-3 ./v_l[:,2],T[:,2],label="PR",linestyle="--",color="g")
plt.semilogx(1e-3 ./v_v[:,2],T[:,2],label="",linestyle="--",color="g")
plt.semilogx(1e-3 ./v_l[:,3],T[:,3],label="PR{TwuAlpha}",linestyle="-.",color="b")
plt.semilogx(1e-3 ./v_v[:,3],T[:,3],label="",linestyle="-.",color="b")
plt.semilogx(1e-3 ./v_l[:,4],T[:,4],label="PR{TwuAlpha, RackettTranslation}",linestyle=(0, (3, 1, 1, 1)),color="purple")
plt.semilogx(1e-3 ./v_v[:,4],T[:,4],label="",linestyle=(0, (3, 1, 1, 1)),color="purple")
plt.semilogx(ρ_l_exp,T_exp,label="Experimental",marker="o",linestyle="",color="k")
plt.semilogx(ρ_v_exp,T_exp,label="",marker="o",linestyle="",color="k")

plt.legend(loc="upper left",frameon=false,fontsize=12) 
plt.xlabel("Density / (mol/L)",fontsize=16)
plt.ylabel("Temperature / K",fontsize=16)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.xlim([1e-4,1e2])
plt.ylim([300,650])
display(plt.gcf())

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Fields

    No fields configured for issues without a type.

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions