This is Part 10 of a series.
Previously, I tested a few flash solvers on mixtures with helium before deciding to make my own reliable solver. In this part I'll look further into whether the predictions of phase equilibria and physical properties can improve by modifying the critical constants of helium used in the GERG-2008 EOS [1]. I will note here that although usually I document things as I try them in this case I have tried some things over the past few weekends and I didn't record them all. I've used the .julia/logs/repl_history.jl to help reconstruct some specifics but mostly I'm relying on memory.
As I now have a more reliable flash solver to work with, I'll begin by (re-)calculating the deviations from the data of the GERG-2008 EOS [1] and two models with modified Tc and ρc from Messerly et al. [2]. The code to do this is based on that from Part 7:
function calculate_statistics(dev)
n = length(dev)
AAD = sum(abs.(dev))/n
Bias = sum(dev)/n
RMS = sqrt(sum(dev.^2)/n)
Max = maximum(dev)>-minimum(dev) ? maximum(dev) : minimum(dev)
return (AAD, Bias, RMS, Max)
end
function flash_statistics(file, N, mixture)
devK1 = zeros(N)
devK2 = zeros(N)
K1 = zeros(N)
K2 = zeros(N)
i = 0
for row in CSV.File(file)
z = [0.5*(row.x1+row.y1), 0.5*(row.x2+row.y2)]
try
(x,y,β) = AdvancedThermo.flash_with_stability(mixture,row.p*1e3,row.T,z)
if β==0.0 || β==1.0
@info "Single phase - $(row.Reference) datum $(row.n)"
else
i += 1
devK1[i] = log(row.y1/row.x1) - log(y[1]/x[1])
devK2[i] = log(row.y2/row.y1) - log(y[2]/x[2])
K1[i] = log(row.y1/row.x1)
K2[i] = log(row.y2/row.x2)
end
catch e
@info "Internal error - $(row.Reference) datum $(row.n)"
println(e.msg)
end
end
(AADK1,BiasK1,RMSK1,MaxK1) = calculate_statistics(devK1[1:i])
(AADK2,BiasK2,RMSK2,MaxK2) = calculate_statistics(devK2[1:i])
println("Statistic lnK(methane) lnK(helium)")
println("AAD $AADK1 $AADK2")
println("Bias $BiasK1 $BiasK2")
println("RMS $RMSK1 $RMSK2")
println("Max $MaxK1 $MaxK2")
return (K1[1:i], K2[1:i], devK1[1:i], devK2[1:i])
end
pTxy_methane_helium.csv
for the various models gives:
Results for GERG-2008 EOS [1]a
-----------------------------
Statistic lnK(methane) lnK(helium)
AAD 0.36743472024844487 2.487039581867481
Bias -0.33229410009211424 -2.4684691120311237
RMS 0.5575148892589821 2.9733531350213824
Max -1.851716348867873 -6.865761787762258
a 3 points predicted as single phase
Results for Messerly A (Tc = 13, ρc = 28400, γv = γT = 1)b
---------------------------------------------------------
Statistic lnK(methane) lnK(helium)
AAD 0.5267775601980111 9.558432793190308
Bias 0.5266834077008744 -9.558432793190308
RMS 0.6184443691443137 10.048008727255144
Max 1.4314091755111868 -20.60754502534418
b 3 points predicted as single phase
Results for Messerly B (Tc = 13, ρc = 28400, γv = 1.086694, γT = 2.044932)
--------------------------------------------------------------------------
Statistic lnK(methane) lnK(helium)
AAD 0.17690224922980394 5.703383457656008
Bias 0.17200694947307213 -5.703383457656008
RMS 0.24685901745331715 6.010484047646062
Max 0.7458249648047901 -13.43577111950453
The results above confirm the conclusions from Part 7,
that using the critical parameters of Messerly et al. [2] with
γv = γT = 1 worsens the predictions
of the methane-helium binary compared to GERG-2008 [1].
The model in which the Messerly et al. critical constants are used
along with interaction parameters based on
arithmetic means
is harder to judge: the prediction of Kmethane improves a bit but
the predictions of Khelium deteriorate substantially. However,
all of the points are correctly classified as comprising two phases which is a small
improvement over GERG.
Even though the prediction of Kmethane values with modified critical parameters is better than GERG-2008, the agreement with data remains noticeably worse than the cubic model with tuned critical constants developed by Rowland et al. (refer Figure 3c of [3]). This suggests to me that I need to fit to the data in order to determine an acceptable set of critical parameters for helium.
As running flash_statistics involves flashing more than 100 points
it takes a while to execute. I attempted to use the flash to find optimal values
of critical temperature and density of classical helium like so:
function fit_critical_constants_flash(file, mixture)
function model(x,p)
mixture.components[2] = new_critical!(mixture.components[2], 1.e6, p[1], p[2]*1e3)
(γv, γT) = arithmetic_γ(mixture.components[1], mixture.components[2])
mixture.γv[1,2] = mixture.γv[2,1] = γv
mixture.γT[1,2] = mixture.γT[2,1] = γT
devK1 = zeros(128)
devK2 = zeros(128)
i = 0
for row in CSV.File(file)
i += 1
z = [0.5*(row.x1+row.y1), 0.5*(row.x2+row.y2)]
try
(x,y,β) = flash_with_stability(mixture,row.p*1e3,row.T,z)
if β==0.0 || β==1.0
# penalise single phase results
devK1[i] = 10.0
devK2[i] = 10.0
else
devK1[i] = log(row.y1/row.x1) - log(y[1]/x[1])
devK2[i] = log(row.y2/row.y1) - log(y[2]/x[2])
end
catch e
# penalise other errors
devK1[i] = 10.0
devK2[i] = 10.0
end
end
return vcat(devK1[1:i], devK2[1:i])
end
p0 = [13.0, 28.0]
fit = LsqFit.curve_fit(model, zeros(256), zeros(256), p0)
return coef(fit)
end
An alternative to running the full flash for every point is to calculate the difference in fugacity between the liquid and vapour phases based on the experimental compositions:
function calculate_goodness_of_fit(Tc,ρc)
methane = AdvancedThermo.get_fluid_component("methane")
helium = AdvancedThermo.get_fluid_component("helium")
# Dummy value for critical pressure
helium = AdvancedThermo.new_critical!(helium,1e6,Tc,ρc)
m = AdvancedThermo.create_mixture([methane,helium])
(γv12, γT12) = AdvancedThermo.arithmetic_γ(methane, helium)
m.γv[1,2] = m.γv[2,1] = γv12
m.γT[1,2] = m.γT[2,1] = γT12
file = "test/pTxy_methane_helium.csv"
ssq = 0.0
for row in CSV.File(file)
T = row.T
p = row.p*1e3
x = [row.x1, row.x2]
lnf_liq = log.(x)+AdvancedThermo.lnϕ(m,AdvancedThermo.ρ(m,p,T,x),T,x)
y = [row.y1, row.y2]
lnf_vap = log.(y)+AdvancedThermo.lnϕ(m,AdvancedThermo.ρ(m,p,T,y),T,y)
#@info lnf_liq, lnf_vap
ssq += sum((lnf_liq.-lnf_vap).^2)
end
rmsd = sqrt(ssq/2/128)
return (ssq, rmsd)
end
arithmetic_γ. An advantage of this method is that it doesn't
need the flash, which actually means I could have tried it prior to
developing the reliable flash. Calculating the goodness-of-fit at a series of
points and plotting the results as a surface can be achieved like so:
using Plots
T = 5:0.1:15.0
ρ = 10000:500:20000
surface(T, ρ, (x,y) -> AdvancedThermo.calculate_goodness_of_fit(x,y)[1])
calculate_goodness_of_fit at several points) which is
a much higher critical temperature than expected. Running the flash_statistics
function gives the following results
Results for Model C (Tc = 101, ρc = 18700, γv = 1.031007, γT = 1.050807)
--------------------------------------------------------------------------
Statistic lnK(methane) lnK(helium)
AAD 0.8794574320713446 3.2290930093654087
Bias -0.7964101645477908 -1.3794488063938029
RMS 2.259885753221556 4.0251045551693805
Max -8.684084051281314 9.607441033426088
Six of the points are predicted as single-phase. Overall, the Kmethane values
are predicted less well compared to the models above with critical parameters from Messerly et al. [2]
while the Khelium values are a bit better predicted. By not performing the
flash, important information about the phase stability is omitted from the goodness-of-fit calculation.
I might dig further into this issue next time.
The code so far can be obtained by downloading the zip file.
The next part of the series is here.
[1] Kunz and Wagner, J. Chem. Eng. Data, 2012, 57, 3032 (link to publisher)
[2] Messerly et al., J. Chem. Eng. Data, 2020, 65, 1028 (link to publisher)
[3] Rowland et al., J. Chem. Eng. Data, 2017, 62, 2799 (link to publisher)
| Go to Home | How to cite this page | Permalink to this page |