Advanced Thermodynamics

Thermodynamic properties of natural gas mixtures containing helium - Part 2

26 June 2021

This piece is Part 2 of a series.

In Part 1, I showed how to calculate the residual contribution to the reduced Helmholtz energy for pure helium. Extending the code to calculate the Helmholtz energy of mixtures is a natural next step.

First I'll add methane as another pure fluid so that there's something to mix with helium.

    # refer Kunz & Wagner (2012) Tables A3 and A5
    np = [
        0.57335704239162
        -0.16760687523730e1
        0.23405291834916
        -0.21947376343441
        0.16369201404128e-1
        0.15004406389280e-1
    ]
    dp = [
        1
        1
        2
        2
        4
        4
    ]
    tp = [
        0.125
        1.125
        0.375
        1.125
        0.625
        1.500
    ]
    ne = [
        0.98990489492918e-1
        0.58382770929055
        -0.74786867560390
        0.30033302857974
        0.20985543806568
        -0.18590151133061e-1
        -0.15782558339049
        0.12716735220791
        -0.32019743894346e-1
        -0.68049729364536e-1
        0.24291412853736e-1
        0.51440451639444e-2
        -0.19084949733532e-1
        0.55229677241291e-2
        -0.44197392976085e-2
        0.40061416708429e-1
        -0.33752085907575e-1
        -0.25127658213357e-2
    ]
    de = [
        1
        1
        1
        2
        3
        6
        2
        3
        3
        4
        4
        2
        3
        4
        5
        6
        6
        7
    ]
    te = [
        0.625
        2.625
        2.750
        2.125
        2.000
        1.750
        4.500
        4.750
        5.000
        4.000
        4.500
        7.500
        14.000
        11.500
        26.000
        28.000
        30.000
        16.000
    ]
    ce = [
        1
        1
        1
        1
        1
        1
        2
        2
        2
        2
        2
        3
        3
        3
        6
        6
        6
        6
    ]
    return FluidComponent(
        name,
        190.564,
        0.0, # NOTE(DR) critical pressure not required (yet)
        10139.342719,
        16.04246e-3,
        np,
        dp,
        tp,
        ne,
        de,
        te,
        ce
    )

For a mixture of fluids, the reduced residual Helmholtz energy α^r is given by Eq (7) of [1]:

The key parts of the mixture model lacking are the composition-dependent reducing functions ρ_r and T_r. The necessary formulae are given by Kunz and Wagner [1] as their Eqs. (16) and (17). These equations are the most complicated to be encountered so far in the GERG-2008 EOS, so some care is needed to implement them. I'm implementing the equations as written in the GERG-2004 Technical Monograph [2] (Eq (7.22a)).

function ρᵣ(x,ρc,βv,γv)
    return inv(vᵣ(x,inv.(ρc),βv,γv))
end
function vᵣ(x,vc,βv,γv)
    result = 0.0
    for j = 1:length(x)
        if x[j]>0.0
            for i = 1:length(x)
                if x[i]>0.0
                    vᵢⱼ = (vc[i]^(1/3)+vc[j]^(1/3))^3/8.0
                    result += ( x[i]*x[j]*βv[i,j]*γv[i,j]*
                                (x[i]+x[j])/(βv[i,j]^2*x[i]+x[j])*vᵢⱼ )
                end
            end
        end
    end
    return result
end
function Tᵣ(x,Tc,βT,γT)
    result = 0.0
    for j = 1:length(x)
        if x[j]>0.0
            for i = 1:length(x)
                if x[i]>0.0
                    Tᵢⱼ = sqrt(Tc[i]*Tc[j])
                    result += ( x[i]*x[j]*βT[i,j]*γT[i,j]*
                                (x[i]+x[j])/(βT[i,j]^2*x[i]+x[j])*Tᵢⱼ )
                end
            end
        end
    end
    return result
end

To test the above, I need the interaction parameters β_v, γ_v, β_T and γ_T for the methane-helium pair. These are given Table A8 of [1]. But first I would like somewhere to store them. For this I'll introduce a new mixture struct and a way to create one from a set of fluids:

mutable struct FluidMixture
    components::Vector{FluidComponent}
    Tc::Vector{Float64}
    ρc::Vector{Float64}
    βv::Matrix{Float64}
    γv::Matrix{Float64}
    βT::Matrix{Float64}
    γT::Matrix{Float64}
end
function create_mixture(fluids::Vector{FluidComponent})
    n = length(fluids)
    Tc = zeros(n)
    ρc = zeros(n)
    for i = 1:n
        Tc[i] = fluids[i].Tc
        ρc[i] = fluids[i].ρc
    end
    (βv, γv, βT, γT) = get_interaction_matrices(fluids)
    return FluidMixture(fluids, Tc, ρc, βv, γv, βT, γT)
end
function get_interaction_matrices(fluids::Vector{FluidComponent})
    n = length(fluids)
    βv = ones(n,n)
    γv = ones(n,n)
    βT = ones(n,n)
    γT = ones(n,n)
    n = length(fluids)
    for j = 1:n
        for i = 1:n
            if fluids[i].Name=="methane" && fluids[j].Name=="helium"
                # only γv and γT non-unity
                γv[i,j] = γv[j,i] = 0.881405683
                γT[i,j] = γT[j,i] = 3.159776855
            else
                # TODO(DR): more interactions as needed
            end
        end
    end
    return (βv, γv, βT, γT)
end

Now I need a function that takes a FluidMixture as input and calculates the residual Helmholtz energy from Eq (7) of [1]. This looks like:

function αʳ(m::FluidMixture,ρ,T,x)
    δ = ρ/ρᵣ(x,m.ρc,m.βv,m.γv)
    τ = Tᵣ(x,m.Tc,m.βT,m.γT)/T
    return sum(x.*αʳ.(δ,τ,m.components))
end
function αʳ(δ,τ,x,m::FluidMixture)
    return sum(x.αʳ(δ,τ,m.components))
end

@testset "αʳ(methane)" begin
    methane = AdvancedThermo.get_fluid_component("methane")
    reltol = 1.e-14
    # benchmark values from TREND 4.0 (GERG-2008)
    αʳ = AdvancedThermo.αʳ(methane,50.0,300.0)
    @test abs(αʳ+2.10257286702117E-03)/αʳ<reltol
    αʳ = AdvancedThermo.αʳ(methane,20000.0,100.0)
    @test abs(αʳ+5.12352321835192)/αʳ<reltol
end
@testset "αʳ(methane+helium)" begin
    methane = AdvancedThermo.get_fluid_component("methane")
    helium = AdvancedThermo.get_fluid_component("helium")
    mixture = AdvancedThermo.create_mixture([methane,helium])
    reltol = 1.e-14
    x = [0.7, 0.3]
    # benchmark values from TREND 4.0 (GERG-2008)
    αʳ = AdvancedThermo.αʳ(mixture,1000.0,300.0,x)
    @test abs(αʳ+1.49341854944043E-02)/αʳ<reltol
    αʳ = AdvancedThermo.αʳ(mixture,20000.0,400.0,x)
    @test abs(αʳ-2.36418364603531E-01)/αʳ<1.e-11
    αʳ = AdvancedThermo.αʳ(mixture,20000.0,200.0,x)
    @test abs(αʳ+5.20295189324906E-01)/αʳ<reltol
    αʳ = AdvancedThermo.αʳ(mixture,20000.0,700.0,x)
    @test abs(αʳ-5.10461374891548E-01)/αʳ<1.e-10
end

I'd call the test of the mixture model above a success, but I'm not quite ready to declare victory. As the mixture parameters β_v and β_T were set to unity for the methane-helium pair, I would like to test another pair for which those parameters differ from 1.0. Consulting [1], nitrogen-helium presents such an opportunity: it's also convenient as I'll need nitrogen before long. As the code to implement the model for nitrogen matches closely that for methane given above I'll not display it here. But it is available in the zip file below, along with new unit tests.

For the nitrogen-helium binary mixture I get the following comparison with TREND 4.0 [3].

@testset "αʳ(nitrogen+helium)" begin
    nitrogen = AdvancedThermo.get_fluid_component("nitrogen")
    helium = AdvancedThermo.get_fluid_component("helium")
    mixture = AdvancedThermo.create_mixture([nitrogen,helium])
    x = [0.7, 0.3]
    # benchmark values from TREND 4.0 (GERG-2008)
    αʳ = AdvancedThermo.αʳ(mixture,1000.0,300.0,x)
    @test abs(αʳ-7.57385600385365E-03)/αʳ<1.e-8
    αʳ = AdvancedThermo.αʳ(mixture,20000.0,400.0,x)
    @test abs(αʳ-4.90263530771741E-01)/αʳ<1.e-12
end

The code so far can be obtained by downloading the zip file. In the next part of the series, I'll develop a way to calculate properties other than the residual Helmholtz energy.

References

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

[2] 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)

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