Advanced Thermodynamics

Reliable flash for mixtures containing helium - Part 1

23 October 2021

In the series about helium-containing mixtures, I created a simple flash procedure for determining vapour-liquid equilibrium in systems such as (methane + helium) and (nitrogen + helium) described by the GERG-2008 Equation of State [1]. I found that the simple solver was not very reliable – sometimes failing to converge or converging to an incorrect solution. I also tested the flash routines from TREND 4.0, REFPROP 10.0 and the julia package ThermoModels.jl and found them all lacking in some respects. In this series of posts, I'd like to arrive at a more reliable method of solving the flash problem for binary mixtures containing helium.

Before writing any code, my first objective is to develop a set of benchmark results. In the course of testing TREND 4.0 [2] and REFPROP 10.0 [3], I flashed (or attempted to flash) the (methane + helium) system at more than 100 distinct pTz conditions with both solvers. Comparing the results, I find that 56 of the calculations are in close agreement (which I chose as being within 2×10⁻⁶ in the computed vapour fraction) which means they probably correspond to the correct equilibrium state. The 56 points are listed in flash_benchmark_GERG2008_methane+helium.csv the first few rows of which look like:

    "T",    "p",     "z1",     "beta",      "y1",        "x1"
    189.01, 6891.2,  0.9086,   0.287420399, 0.784028555, 0.95884614
    189.01, 4974.5,  0.97685,  0.295109999, 0.950485776, 0.987887674
    188.01, 4819.4,  0.9757,   0.32786611,  0.94800237,  0.989210871

The following code runs the flash calculations and compares to the benchmark results:

@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 = "flash_benchmark_GERG2008_methane+helium.csv"
    correct = 0
    errors = 0
    for row in CSV.File(file)
        z = [row.z1, 1-row.z1]
        try 
            # Uncomment one line below
            #(x,y,β) = AdvancedThermo.flash(mixture,row.p*1e3,row.T,z)
            (x,y,β) = AdvancedThermo.flash_helium_binary(mixture,row.p*1e3,row.T,z)
            #(x,y,β) = AdvancedThermo.flash_initial_K(mixture,row.p*1e3,row.T,z,[row.y1/row.x1, (1-row.y1)/(1-row.x1)])
            if β==0.0 || β==1.0
                @info "Single phase - row $(row.n)"
            elseif abs(β-row.beta)>2.e-6
                @info "Incorrect β - row $(row.n)"
            elseif abs(x[1]-row.x1)>2.e-6
                @info "Incorrect x1 - row $(row.n)"
            elseif abs(y[1]-row.y1)>2.e-6
                @info "Incorrect y1 - row $(row.n)"
            else
                correct += 1
            end
        catch e
            @info "Internal error - row $(row.n)"
            println(e.msg)
            errors += 1
        end
    end
    @info "Correct: $correct,  Incorrect: $(56-correct-errors),  Errors: $errors"
end

    flash()                    Correct: 2,   Incorrect: 54,  Errors: 0
    flash_helium_binary()      Correct: 37,  Incorrect: 19,  Errors: 0
    flash_initial_K()          Correct: 46,  Incorrect: 10,  Errors: 0

I'll improve the flash using some of the techniques described by Michelsen and Mollerup [4]. First I'll implement a stability analysis based on finding minima in the tangent plane distance. The code looks like the following:

function stability_helium_binary(m::FluidMixture,p,T,z)
    lnK = lnK_helium_binary(m,p,T)
    K = exp.(lnK)
    W = K.*z
    Wv = minimise_tm(m,p,T,z,W)
    W = z./K
    Wl = minimise_tm(m,p,T,z,W)
    tpd_v = tpd(m,p,T,z,Wv./sum(Wv))
    tpd_l = tpd(m,p,T,z,Wl./sum(Wl))
    tm_v = tm(m,p,T,z,Wv)
    tm_l = tm(m,p,T,z,Wl)
    if tpd_v < tpd_l
        return (Wv, tpd_v, tm_v)
    else
        return (Wl, tpd_l, tm_l)
    end
end
function minimise_tm(m::FluidMixture,p,T,z,W)
    d = log.(z) .+ lnϕ(m,ρ(m,p,T,z),T,z)
    for i = 1:200
        w = W./sum(W)
        # Michelsen/Mollerup (10.51)
        Wnew = exp.(d .- lnϕ(m,ρ(m,p,T,w),T,w))
        if maximum(abs.(Wnew.-W))<1.e-5
            return Wnew
        end
        W .= Wnew
    end
    @warn "failed to minimise tm(W)"
    return W
end
function tpd(m::FluidMixture,p,T,z,w)
    # Michelsen/Mollerup (9.4) and (10.42)
    return sum(w.*(log.(w).+lnϕ(m,ρ(m,p,T,w),T,w).-log.(z).-lnϕ(m,ρ(m,p,T,z),T,z)))
end
function tm(m::FluidMixture,p,T,z,W)
    WT = sum(W)
    w = W./WT
    return 1 - WT + WT*log(WT) + WT*tpd(m,p,T,z,w)
end

@testset "stability 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]
        p = row.p*1.e3
        T = row.T
        try
            (W,tpd,tm) = AdvancedThermo.stability_helium_binary(mixture,p,T,z)
            if tm >= 0.0
                @info "tm>0 - row $(row.n)"
            end
        catch e
            @info "Internal error - row $(row.n)"
        end
    end
end

    "n", "T",    "p",     "z1",     "beta",      "y1",        "x1"
    32,  154.01, 5560.6,  0.63655,  0.515872068, 0.309908188, 0.984609627

    minimise_tm(m, 5560.6e3, 154.01, [0.63655, 1-0.63655], 
                        [2.889229990105307, 1.5737191000899076e-5])

#  result of "@info Wnew, tmW"
... snip first few iterations ...
[ Info: ([1.4812301176977118, 0.009811385532533598], -2.5764174126546933)
[ Info: ([8.710330873755286, 1.3834227689297245e-26], 7.721212493767755)
[ Info: ([1.4812301176977118, 0.009811385532533598], -2.5764174126546933)
[ Info: ([8.710330873755286, 1.3834227689297245e-26], 7.721212493767755)

For stability calculations the convergence behaviour is dictated by eigenvalues of the matrix containing partial derivatives like [5]:

• n_i(∂lnφ_i /∂n_j)

function ∂lnϕ∂n(m::FluidMixture,V,T,n)
    dna = AdvancedThermo.derivatives_of_nαʳ(m,V,T,n)
    dp = AdvancedThermo.derivatives_of_p(m,V,T,n)
    ∂p∂V = DiffResults.gradient(dp)[1]
    ∂p∂n = DiffResults.gradient(dp)[3:end]
    ∂lnϕ∂n = DiffResults.hessian(dna)[3:end,3:end] .+ 1 .+
                  ∂p∂n*∂p∂n'./(R*T*∂p∂V)
    return ∂lnϕ∂n
end

I examined parts of the Gibbs energy curve and tangent plane distance with code like

    for w=0.01:0.01:0.99
        ρ = AdvancedThermo.ρ(mixture, 5560.6e3, 154.01, [w, 1-w])
        println(AdvancedThermo.g(mixture, ρ, 154.01, [w, 1-w]))
    end
    for w=0.01:0.01:0.99
        println(AdvancedThermo.tpd(mixture, 5560.6e3, 154.01, [0.63655, 1-0.63655], [w, 1-w]))
    end

function minimise_tm(m::FluidMixture,p,T,z,W)
    d = log.(z) .+ lnϕ(m,ρ(m,p,T,z),T,z)
    tm_ = tm(m,p,T,z,W)
    for i = 1:200
        w = W./sum(W)
        # Michelsen/Mollerup (10.51)
        Wnew = exp.(d .- lnϕ(m,ρ(m,p,T,w),T,w))
        tmnew = tm(m,p,T,z,Wnew)
        if maximum(abs.(Wnew.-W))<1.e-5
            return Wnew
        end
        # reject steps that don't decrease tm
        if tmnew >= tm_
            return W
        end
        W .= Wnew
        tm_ = tmnew
    end
    @warn "failed to minimise tm(W)"
    return W
end

The code so far can be obtained by downloading the zip file.

References

[1] Kunz and Wagner, J. Chem. Eng. Data, 2012, 57, 3032 (link to publisher)

[2] Span et al., TREND. Thermodynamic Reference and Engineering Data 4.0. Lehrstuhl für Thermodynamik, Ruhr-Universität Bochum, 2019

[3] Lemmon et al., REFPROP Documentation, Release 10.0, NIST (link to FAQ)

[4] Michelsen and Mollerup, Thermodynamic Models: Fundamentals & Computational Aspects, 2nd Edition, Tie-Line Publications, 2007

[5] Heidemann and Michelsen, Ind. Eng. Chem. Res., 1995, 34, 958 (link to publisher)