Advanced Thermodynamics

Thermodynamic properties of natural gas mixtures containing helium - Part 5

7 July 2021

This is Part 5 of a series.

In Part 4, I implemented a density solver and function for calculating the fugacity coefficients in a mixture. Both of these will be used heavily in this part, which focuses on implementing a phase equilibrium solver.

In the classic 2-phase equilibrium problem, given pressure p, temperature T and overall/feed composition z, the task is to find a liquid composition x, vapour composition y, and vapour phase fraction β such that for all components i:

• (1-β)x_i + βy_i = z_i, and
• x_iφ_i,liq = y_iφ_i,vap

1. Estimate K_i (≡ y_i/x_i)
2. Solve the Rachford-Rice equation for β, then calculate x_i and y_i
3. Calculate K_i = φ_i,liq(p,T,x)/φ_i,vap(p,T,y)
4. Repeat steps 2 and 3 until convergence

For step 1 above, the most common general method for calculating K_i is due to Wilson. The function can be implemented like so:

function lnK_Wilson(m::FluidMixture,p,T)
    n = length(m.components)
    lnK = zeros(n)
    for i = 1:n
        lnK[i] = log(m.components[i].pc/p)+5.37*(1+m.components[i].ω)*(1-m.components[i].Tc/T)
    end
    return lnK
end

For step 2, the Rachford-Rice equation and it's solver can be written succinctly:

function solve_xyβ(z,K)
    (f0,df0) = RachfordRice(z,K,0.0)
    (f1,df1) = RachfordRice(z,K,1.0)
    if f0 > 0.0 && f1 < 0.0
        β = 0.5
        (f,df) = RachfordRice(z,K,β)
        while abs(f)>1.e-5
            β -= f/df
            (f,df) = RachfordRice(z,K,β)
        end
        x = @.(z/(1+β*(K-1)))
        y = K.*x
        return (x,y,β)
    else
        return (z,z, f0 > 0.0 ? 1.0 : 0.0)
    end
end
function RachfordRice(z,K,β)
    s = @.((K-1)/(1+β*(K-1)))
    return (z's, -z's.^2)
end

@testset "Rachford-Rice solver" begin
    # adapted from 'Advanced Thermodynamics - Lecture 16 (2019)'
    z = [0.0191,0.0433,0.7046,0.1851,0.0479]
    K = [0.434,8.53,2.9696,0.3038,0.0593]
    (x,y,β) = AdvancedThermo.solve_xyβ(z,K)
    x_ref = [0.0368, 0.0058, 0.2633, 0.4541, 0.2399]
    y_ref = [0.0160, 0.0499, 0.7820, 0.1379, 0.0142]
    β_ref = 0.85078
    tol = 1.e-4
    for i = 1:length(x_ref)
        @test abs(x[i]-x_ref[i])<tol
        @test abs(y[i]-y_ref[i])<tol
    end
    @test abs(β-β_ref)<tol
end

For step 3, last time I developed a function for calculating φ_i(ρ,T,x) and a density solver. Using both of these and the code above, the full algorithm looks like:

function flash(m::FluidMixture,p,T,z)
    lnK = lnK_Wilson(m,p,T)
    K = exp.(lnK)
    β_prev = -1.0
    (x,y,β) = solve_xyβ(z,K)
    while abs(β-β_prev)>1.e-6
        lnϕ_liq = lnϕ(m,ρ(m,p,T,x),T,x)
        lnϕ_vap = lnϕ(m,ρ(m,p,T,y),T,y)
        β_prev = β
        (x,y,β) = solve_xyβ(z,exp.(lnϕ_liq.-lnϕ_vap))
    end
    return (x,y,β)
end

    @testset "pT flash (methane+helium)" begin
    methane = AdvancedThermo.get_fluid_component("methane")
    helium = AdvancedThermo.get_fluid_component("helium")
    mixture = AdvancedThermo.create_mixture([methane,helium])
    (x,y,β) = AdvancedThermo.flash(mixture,2.e6,90.0,[0.5,0.5])
    # benchmark values from TREND 4.0
    x_ref = [0.999896, 0.000104]
    y_ref = [0.007754, 0.992246]
    β_ref = 0.503856
    tol = 1.e-6
    for i = 1:length(x_ref)
        @test abs(x[i]-x_ref[i])<tol
        @test abs(y[i]-y_ref[i])<tol
    end
    @test abs(β-β_ref)<tol
end

@testset "pT flash (methane+helium)" begin
    methane = AdvancedThermo.get_fluid_component("methane")
    helium = AdvancedThermo.get_fluid_component("helium")
    mixture = AdvancedThermo.create_mixture([methane,helium])
    (x,y,β) = AdvancedThermo.flash(mixture,3.e6,90.0,[0.5,0.5])
    # benchmark values from TREND 4.0
    x_ref = [0.999847, 0.000153]
    y_ref = [0.006025, 0.993975]
    β_ref = 0.502955
    tol = 1.e-6
    for i = 1:length(x_ref)
        @test_broken abs(x[i]-x_ref[i])<tol
        @test_broken abs(y[i]-y_ref[i])<tol
    end
    @test_broken abs(β-β_ref)<tol
end

Looking at the failure case, the final K values are around 0.006/0.9998 = 0.006 for methane and 0.994/0.00015 = 6600 for helium. In contrast, the K values from Wilson's method at this condition are 0.0035 and 1.61 respectively. The K value for helium is particularly poor, and this may help explain the failure to flash the mixture correctly. As the equilibrium phases are nearly pure, better estimates of the K values can be made by calculating the partial fugacity coefficients assuming a perfect split of the components. This is incorporated into a flash algorithm below:

function flash_helium_binary(m::FluidMixture,p,T,z)
    x1 = [1.0, 0.0]
    x2 = [0.0, 1.0]
    ρ1 = ρ(m,p,T,x1)
    ρ2 = ρ(m,p,T,x2)
    lnϕ_1 = lnϕ(m,ρ1,T,x1)
    lnϕ_2 = lnϕ(m,ρ2,T,x2)
    K = ρ1>ρ2 ? exp.(lnϕ_1.-lnϕ_2) : exp.(lnϕ_2.-lnϕ_1)
    return flash_initial_K(m,p,T,z,K)
end
function flash_initial_K(m::FluidMixture,p,T,z,K)
    β_prev = -1.0
    (x,y,β) = solve_xyβ(z,K)
    nit = 0
    while abs(β-β_prev)>1.e-6 && nit<100
        nit = nit + 1
        lnϕ_liq = lnϕ(m,ρ(m,p,T,x),T,x)
        lnϕ_vap = lnϕ(m,ρ(m,p,T,y),T,y)
        β_prev = β
        (x,y,β) = solve_xyβ(z,exp.(lnϕ_liq.-lnϕ_vap))
    end
    return (x,y,β)
end

To more thoroughly test the flash, I'll use the same 128 data points for the (methane + helium) binary system that were used in Rowland et al. [3]. These data originate from Heck and Hiza [4], DeVaney et al. [5], and Rhodes et al. [6]. The data are included in the zip file (see below) and can be read with CSV.jl then sent to the flash:

@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_initial_K(mixture,row.p*1e3,row.T,z,K)
            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

    [ Info: Internal error - Heck and Hiza, 1969 datum 4
    log will only return a complex result if called with a complex argument. Try log(Complex(x)).
    [ Info: Single phase - Heck and Hiza, 1969 datum 5
    [ Info: Single phase - Heck and Hiza, 1969 datum 11
    [ Info: Single phase - Heck and Hiza, 1969 datum 13
    [ Info: Single phase - Heck and Hiza, 1969 datum 21
    [ Info: Single phase - Heck and Hiza, 1969 datum 28
    [ Info: Single phase - Heck and Hiza, 1969 datum 30
    [ Info: Internal error - DeVaney et al., 1971 datum 4
    log will only return a complex result if called with a complex argument. Try log(Complex(x)).
    [ Info: Internal error - DeVaney et al., 1971 datum 5
    log will only return a complex result if called with a complex argument. Try log(Complex(x)).
    [ Info: Single phase - Rhodes et al., 1971 datum 5
    [ Info: Single phase - Rhodes et al., 1971 datum 35
    [ Info: Single phase - Rhodes et al., 1971 datum 44

It seems there are 2 options for what to tackle next:

1. debug the errors in the flash and get the number of failures to zero, or
2. try running the pTxy test cases with modified critical parameters

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

References

[1] Kunz et al., GERG Technical Monograph 15: The GERG-2004 Wide-Range Equation of State for Natural Gases and Other Mixtures, Groupe Européen de Recherches Gazières, 2007 (link to PDF)

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

[3] Rowland et al., J. Chem. Eng. Data, 2017, 62, 2799 (link to publisher)

[4] Heck and Hiza, AIChE J., 1967, 13, 593 (link to publisher)

[5] DeVaney et al., J. Chem. Eng. Data, 1971, 16, 158 (link to publisher)

[6] Rhodes et al., J. Chem. Eng. Data, 1971, 16, 19 (link to publisher)