In Part 2, I improved the density solver which led to passing most of the flash benchmarks. However, there were a few cases where my flash algorithm made different predictions compared to the benchmarks so I want to look more closely at these points.
The largest deviation occurred for point 2 of the benchmark — around 0.0004 in terms of vapour fraction. The benchmark for this point looks like
"n", "T", "p", "z1", "beta", "y1", "x1"
2 , 189.01, 4974.5, 0.97685, 0.295109999, 0.950485776, 0.987887674
which can be flashed by running
(x,y,β) = AdvancedThermo.flash_with_stability(m,4974.5e3,189.01,[0.97685, 1-0.97685])
([0.9878799466525535, 0.012118742819828115], [0.9504574971507432, 0.04954563868526698], 0.2947414399482296)I think the first thing to check is whether the fugacities are equal in the vapour and liquid. The fugacity (as natural log) can be obtained for the liquid via
lnf_liq = log.(x) .+ AdvancedThermo.lnϕ(m,AdvancedThermo.ρ(m,p,T,x),T,x)
| Component | lnfliq | lnfvap |
|---|---|---|
| methane | -0.505120231268394 | -0.5051202313362685 |
| helium | -1.9822780366836952 | -1.9822780354073914 |
| Component | lnfliq | lnfvap |
|---|---|---|
| methane | -0.50511130782736074 | -0.5051113079537916 |
| helium | -1.9826009856476388 | -1.982600983002531 |
g = ( β*AdvancedThermo.g(m,AdvancedThermo.ρ(m,p,T,y),T,y)
+ (1-β)*AdvancedThermo.g(m,AdvancedThermo.ρ(m,p,T,x),T,x) )
I suspect the most likely problem is that the numerical tolerances in some of my tests are inadequate. There are only a few places where tolerances are needed, these are:
abs(f)<1.e-5 has been sufficient for termination.
Changing this to abs(f)<1.e-12 gives
the following result for the flash of benchmark point 2:
x = [0.9878876739951483, 0.012112326004851614] y = [0.9504857761550682, 0.049514223844931746] β = 0.29510999795631604which is enough of a shift to produce excellent agreement with the benchmark result :). Running through the entire benchmark now gives the summary
[ Info: Correct: 56, Incorrect: 0, Errors: 0i.e. all of the benchmark cases are passing successfully! To help ensure it stays that way, I've added a set of unit tests based on flashing the benchmark cases.
@testset "flash benchmark (methane+helium)" begin
methane = AdvancedThermo.get_fluid_component("methane")
helium = AdvancedThermo.get_fluid_component("helium")
mixture = AdvancedThermo.create_mixture([methane,helium])
file = "flash_benchmark_GERG2008_methane+helium.csv"
for row in CSV.File(file)
z = [row.z1, 1-row.z1]
(x,y,β) = AdvancedThermo.flash_with_stability(mixture,row.p*1e3,row.T,z)
@test abs(β-row.beta)<2.e-6
@test abs(x[1]-row.x1)<2.e-6
@test abs(y[1]-row.y1)<2.e-6
end
end
I also checked one of the points from the full data set that gave a single phase result last time ("Rhodes et al., 1971 datum 35") and it does seem that stability analysis is converging slowly to the trivial solution. That, combined with the results from other software and the figures in this paper [1], convinces me that a single phase is the correct flash outcome for that point with the GERG-2008 EOS.
I won't be shocked to re-visit the flash algorithm in the future, but for now I believe the flash for helium mixtures I've developed surpasses all of the available flashes I tried, so I think it's time to put this series on hold and look once again at the problem of improving the prediction of properties for mixtures containing helium.
The code so far — including the flash benchmarks — can be obtained by downloading the zip file.
[1] 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 |