In Part 1, I prepared some benchmark results for flash calculations on mixtures comprising (methane + helium) and developed a stability analysis function for those mixtures. In this part I will continue developing the flash algorithm and hopefully get to a point where it passes all of the benchmarks.
My simple idea for a flash algorithm is to use stability analysis to initialise the K values and then proceed with successive substitution. This looks like:
function flash_with_stability(m::FluidMixture,p,T,z)
(W, K, tpd, tm) = stability_helium_binary(m,p,T,z)
if tpd < -1.e-12
(x,y,β) = flash_initial_K(m,p,T,z,K)
else
x = z
y = z
β = 0.0
end
return (x,y,β)
end
[ Info: Incorrect β - row 1: 0.2875542865072627 vs. 0.287420399
[ Info: Incorrect β - row 2: 0.294740164238877 vs. 0.295109999
[ Info: Incorrect β - row 3: 0.3277482268417944 vs. 0.32786611
[ Info: Incorrect β - row 4: 0.3896284941313747 vs. 0.389634285
[ Info: Incorrect β - row 6: 0.3272473661420213 vs. 0.441713204
[ Info: Incorrect β - row 10: 0.47201090595788686 vs. 0.472021861
[ Info: Single phase - row 34
[ Info: Incorrect β - row 45: 0.5420315210175268 vs. 0.542033639
[ Info: Incorrect β - row 48: 0.47235894526563854 vs. 0.565106308
[ Info: Incorrect β - row 49: 0.488235844790858 vs. 0.571613513
[ Info: Incorrect β - row 50: 0.5883414886923838 vs. 0.588345059
[ Info: Incorrect β - row 52: 0.5994749936780771 vs. 0.599477961
[ Info: Incorrect β - row 53: 0.6004169917750538 vs. 0.600420662
[ Info: Incorrect β - row 55: 0.6358670974413189 vs. 0.635869432
[ Info: Correct: 42, Incorrect: 14, Errors: 0
As discussed previously,
the tolerance in the benchmark is set to 2×10−6
in the vapour fraction or phase composition. Most of the results (42 out of 56)
are within this level of agreement. Although this leaves 14 incorrect results,
the results show that many of the incorrect results are only slightly outside
the tolerance so there may be simple explanations for the deviation.
However, there are some larger deviations (rows 6, 48 and 49) and one
point is wrongly reported as single phase. I will start debugging by looking
at this point.
Code to run the flash for the point in question looks like:
AdvancedThermo.flash_with_stability(m,6444.3e3,124.85,[0.521855, 1-0.521855])
@show directive placed in
flash_initial_K() I get the following information:
(x, y, β) = solve_xyβ(z, K) = ([0.9918795159975649, 0.008114755153459068],
[0.19613033246216136, 0.803873637601877], 0.5906691778297667)
(x, y, β) = solve_xyβ(z, exp.(lnϕ_liq .- lnϕ_vap)) = ([0.9943226251521927, 0.005677367740930493],
[0.09078587196911422, 0.9092141345150448], 0.522909138435965)
(x, y, β) = solve_xyβ(z, exp.(lnϕ_liq .- lnϕ_vap)) = ([0.521855, 0.47814500000000004],
[0.521855, 0.47814500000000004], 1.0)
(x, y, β) = solve_xyβ(z, exp.(lnϕ_liq .- lnϕ_vap)) = ([0.521855, 0.47814500000000004],
[0.521855, 0.47814500000000004], 1.0)
I can see that after 1 iteration of successive substitution the vapour fraction
is already quite close to the expected solution (0.523 vs 0.518) but on the
next iteration the phase fraction jumps to 1.0. To find out the cause for the jump
I check the values of lnϕ_vap and lnϕ_liq and find the
following:
lnϕ_liq = lnϕ(m, ρ(m, p, T, x), T, x) = [-2042.6099842016192, 52630.23469784307]
lnϕ_vap = lnϕ(m, ρ(m, p, T, y), T, y) = [-0.5370678892404568, 0.0861504737620466]
These values of lnϕ_liq are completely unrealistic. I suspect that
the density solver is to blame (and is probably responsible for the convergence
issue I was looking at in Part 1).
The density that is being reported from
AdvancedThermo.ρ(m, p, T, [0.9943226251521927, 0.005677367740930493])
To help the solver avoid the non-physical region, I'll stop the search if it looks like the next step will pass a local minimum/maximum in density. The search for a root in the vapour region now looks like:
ρ_lo = 1.e-6
(p_lo,∂p∂ρ,∂²p∂ρ²) = derivatives_of_p_wrt_ρ(m,ρ_lo,T,x)
nit = 0
while abs(p_lo-p)/p > 1.e-10 && nit < 30
dρ = (p-p_lo)/∂p∂ρ
dρ_maximum = -∂p∂ρ/∂²p∂ρ² > 0 ? -∂p∂ρ/∂²p∂ρ² : 1e100
if ∂p∂ρ > 0.0 && dρ < dρ_maximum
if ρ_lo + dρ > 0.0
ρ_lo += dρ
(p_lo,∂p∂ρ,∂²p∂ρ²) = derivatives_of_p_wrt_ρ(m,ρ_lo,T,x)
else
ρ_lo = -1.0
p_lo = p
end
else
ρ_lo = -1.0
p_lo = p
end
nit += 1
end
ρ_hi = 1.e6
(p_hi,∂p∂ρ,∂²p∂ρ²) = derivatives_of_p_wrt_ρ(m,ρ_hi,T,x)
nit = 0
while abs(p_hi-p)/p > 1.e-10 && nit < 30
dρ = (p-p_hi)/∂p∂ρ
dρ_minimum = -∂p∂ρ/∂²p∂ρ²
if dρ>dρ_minimum
ρ_hi += dρ
(p_hi,∂p∂ρ,∂²p∂ρ²) = derivatives_of_p_wrt_ρ(m,ρ_hi,T,x)
else
ρ_hi = -1.0
p_hi = p
end
nit += 1
end
@testset "density solver (methane+helium)" begin
methane = AdvancedThermo.get_fluid_component("methane")
helium = AdvancedThermo.get_fluid_component("helium")
mixture = AdvancedThermo.create_mixture([methane,helium])
z = [0.9943226251521927, 0.005677367740930493]
@test abs(AdvancedThermo.ρ(mixture, 6.4443e6, 124.85, z)-25533.416)<0.001
@test abs(AdvancedThermo.ρ(mixture, 6.4443e6, 124.85, [0.090786, 0.909214])-5985.41)<0.01
@test abs(AdvancedThermo.ρ(mixture, 2e8, 125.0, [0.99, 0.01])-31357.3)<0.1
@test abs(AdvancedThermo.ρ(mixture, 1e9, 125.0, [0.99, 0.01])-38646.3)<0.1
end
julia> AdvancedThermo.flash_with_stability(m,6444.3e3,124.85,[0.521855, 1-0.521855])
([0.9945434613197645, 0.0054565362636050645], [0.08241949796700944, 0.9175805042796146], 0.5182283114043752)
This vapour fraction agrees with the benchmark value to within around 2×10−8
which is more than adequate. The summary from flashing all of the
benchmark cases is now looking quite good too:
[ Info: Incorrect β - row 1: 0.2875542865053227 vs. 0.287420399
[ Info: Incorrect β - row 2: 0.29474016411519666 vs. 0.295109999
[ Info: Incorrect β - row 3: 0.3277482268417953 vs. 0.32786611
[ Info: Incorrect β - row 4: 0.38962849307285896 vs. 0.389634285
[ Info: Incorrect β - row 10: 0.47201090595788686 vs. 0.472021861
[ Info: Incorrect β - row 45: 0.5420315210175268 vs. 0.542033639
[ Info: Incorrect x1 - row 48: 0.9577496218023707 vs. 0.957751919
[ Info: Incorrect β - row 49: 0.5716106842186449 vs. 0.571613513
[ Info: Incorrect β - row 50: 0.5883414886883639 vs. 0.588345059
[ Info: Incorrect β - row 52: 0.5994749936780758 vs. 0.599477961
[ Info: Incorrect β - row 53: 0.6004169917750543 vs. 0.600420662
[ Info: Incorrect β - row 55: 0.6358670974413181 vs. 0.635869432
[ Info: Correct: 44, Incorrect: 12, Errors: 0
Most of the "incorrect" cases look to be falling only just outside
the chosen tolerance. However, a couple of points deviate by much more
than the tolerance which warrants some investigation. I also checked the
flash against the full data set by running:
@testset "pTxy data (methane+helium)" begin
methane = AdvancedThermo.get_fluid_component("methane")
helium = AdvancedThermo.get_fluid_component("helium")
mixture = AdvancedThermo.create_mixture([methane,helium])
file = "pTxy_methane_helium.csv"
for row in CSV.File(file)
z = [0.5*(row.x1+row.y1), 0.5*(row.x2+row.y2)]
K = [row.y1/row.x1, row.y2/row.x2]
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)"
end
catch e
@info "Internal error - $(row.Reference) datum $(row.n)"
println(e.msg)
end
end
end
┌ Warning: failed to minimise tm(W) └ @ AdvancedThermo C:\Jess88\JesWeb\thermo\dev\AdvancedThermo\src\AdvancedThermo.jl:75 [ Info: Single phase - Rhodes et al., 1971 datum 35 ┌ Warning: failed to minimise tm(W) └ @ AdvancedThermo C:\Jess88\JesWeb\thermo\dev\AdvancedThermo\src\AdvancedThermo.jl:75 ┌ Warning: failed to minimise tm(W) └ @ AdvancedThermo C:\Jess88\JesWeb\thermo\dev\AdvancedThermo\src\AdvancedThermo.jl:75 [ Info: Single phase - Rhodes et al., 1971 datum 43 [ Info: Single phase - Rhodes et al., 1971 datum 44so these three points also deserve some attention. I will note that both TREND [1] and REFPROP [2] appear to consider these points as single phase when modelled using the GERG-2008 EOS [3] so it might not indicate a problem with the flash algorithm.
The code so far can be obtained by downloading the zip file. In the next part I'll examine more closely the disagreements with the flash benchmarks.
[1] Span et al., TREND. Thermodynamic Reference and Engineering Data 4.0. Lehrstuhl für Thermodynamik, Ruhr-Universität Bochum, 2019
[2] Lemmon et al., REFPROP Documentation, Release 10.0, NIST (link to FAQ)
[3] Kunz and Wagner, J. Chem. Eng. Data, 2012, 57, 3032 (link to publisher)
| Go to Home | How to cite this page | Permalink to this page |