Advanced Thermodynamics

Thermodynamic properties of natural gas mixtures containing helium - Part 6

9 July 2021

This is Part 6 of a series.

In Part 5, I implemented a basic flash algorithm and applied it to around 130 pTxy literature data for (methane + helium), with 90% of cases giving acceptable results. Before addressing the failures in the flash, I want to investigate modifying the critical constants of helium and how that affects the Helmholtz EOS.

In [1], we found that changing the critical constants of helium in the Peng-Robinson EOS improved significantly the predictions of fluid phase equilibria of mixtures involving helium. In [2], Messerly et al. calculated the critical parameters of 'classical' helium but didn't investigate the effect on mixture predictions. The alternate parameter sets are given below.

Critical density and critical temperature are key quantities in the GERG-2008 EOS [3]. In the first part of this series I implemented models for the Helmholtz energy of pure fluids from the GERG-2008 EOS. For pure fluids, these models depend on critical parameters through δ = ρ/ρ_c and τ = T_c/T. The k-th polynomial-type term contributes to the Helmholtz energy like this:

    n_k × (ρ/ρ_c)^d_k × (T_c/T)^t_k

If a different set of critical parameters are used, say ρ_e and T_e, new EOS parameters (n_k, d_k, t_k) must be found. However, altering d_k or t_k would alter the functional form and mean that the amounts contributed to the Helmholtz energy would differ at the same ρ and T. Therefore the only change must come from n_k. I find the new coefficients n*_k by rearranging the following equation

    n_k × (ρ/ρ_c)^d_k × (T_c/T)^t_k = n*_k × (ρ/ρ_e)^d_k × (T_e/T)^t_k

to arrive at a relationship giving the new coefficients in terms of old:

    n*_k = n_k × (ρ_e/ρ_c)^d_k × (T_c/T_e)^t_k

The other key expressions in the residual Helmholtz energy are of exponential-type:

    n_k × (ρ/ρ_c)^d_k × (T_c/T)^t_k × exp[−(ρ/ρ_c)^c_k]

As before, the introduction of new critical parameters must not change the Helmholtz energy. Due to the presence of the exponential function, this is more complicated than above and requires introducing a new parameter γ*_k:

    n*_k × (ρ/ρ_e)^d_k × (T_e/T)^t_k × exp[−γ*_k(ρ/ρ_e)^c_k]

where n*_k and n_k are related as above and γ*_k = (ρ_e/ρ_c)^c_k.

To accomodate this new parameter, I'll add it to the FluidComponent struct and set it to unity for unmodified fluids. Then the function for the reduced residual Helmholtz energy due to exponential-type terms changes from

function αʳ(δ,τ,n,d,t,c)
    return sum(@.(n*δ^d*τ^t*exp(-δ^c)))
end

function αʳ(δ,τ,n,d,t,γ,c)
    return sum(@.(n*δ^d*τ^t*exp(-γ*δ^c)))
end

The final pure component parameters to deal with occur in the ideal contribution to the Helmholtz energy. For helium, the functional form is simple [3]:

log(ρ/ρ_c) + n₁ + n₂T_c/T + n₃log(T_c/T)

which will become

log(ρ/ρ_e) + n*₁ + n*₂T_e/T + n₃log(T_e/T)

upon substitution of n*₁ = log(ρ_e/ρ_c) + n₁ + n₃log(T_c/T_e), and n*₂ = n₂T_c/T_e.

Fluids other than helium have a more complicated ideal contribution to the Helmholtz energy, including terms like [3]

    n_k × ln[|sinh(v_kT_c/T)|]

These are dealt with by setting v*_k = v_kT_c/T_e and leaving n_k as-is. At last, the function that sets new critical parameters and adjusts the parameters in the EOS for a pure fluid looks like:

function new_critical!(fluid::FluidComponent,pc,Tc,ρc)
    ρr = ρc/fluid.ρc
    Tr = fluid.Tc/Tc
    fluid.pc = pc
    fluid.Tc = Tc
    fluid.ρc = ρc
    # ideal contribution
    fluid.ni[1] += log(ρr) + fluid.ni[3]*log(Tr)
    fluid.ni[2] *= Tr
    @. fluid.vi *= Tr
    # residual contribution
    @. fluid.np *= ρr^fluid.dp * Tr^fluid.tp
    @. fluid.ne *= ρr^fluid.de * Tr^fluid.te
    @. fluid.γe *= ρr^fluid.ce
    return fluid
end

@testset "critical parameters (helium)" begin
    T = 300.0
    p = 1.e5
    x = [1.0]
    helium = AdvancedThermo.get_fluid_component("helium")
    m = AdvancedThermo.create_mixture([helium])
    da = AdvancedThermo.derivatives_of_a(m,AdvancedThermo.ρ(m,p,T,x),T,x)
    
    pc = 0.93e6
    Tc = 13.0
    ρc = 28.4e3
    helium_new = AdvancedThermo.get_fluid_component("helium")
    helium_new = AdvancedThermo.new_critical!(helium_new,pc,Tc,ρc)
    m_new = AdvancedThermo.create_mixture([helium_new])
    da_new = AdvancedThermo.derivatives_of_a(m_new,AdvancedThermo.ρ(m_new,p,T,x),T,x)
    
    # compare Helmholtz energy and derivatives (new vs. old)
    tol = 1.e-13
    val = DiffResults.value(da)
    val_new = DiffResults.value(da_new)
    grad = DiffResults.gradient(da)
    grad_new = DiffResults.gradient(da_new)
    hess = DiffResults.hessian(da)
    hess_new = DiffResults.hessian(da_new)
    @test abs((val-val_new)/val)<tol
    @test abs((grad[1]-grad_new[1])/grad[1])<tol
    @test abs((grad[2]-grad_new[2])/grad[2])<tol
    @test abs((grad[3]-grad_new[3])/grad[3])<tol
    @test abs((hess[1,1]-hess_new[1,1])/hess[1,1])<tol
    @test abs((hess[1,2]-hess_new[1,2])/hess[1,2])<tol
    @test abs((hess[1,3]-hess_new[1,3])/hess[1,3])<tol
    @test abs((hess[2,2]-hess_new[2,2])/hess[2,2])<tol
    @test abs((hess[2,3]-hess_new[2,3])/hess[2,3])<tol
    @test abs((hess[3,3]-hess_new[3,3])/hess[3,3])<tol
end
@testset "critical parameters (methane)" begin
    T = 300.0
    p = 1.e5
    x = [1.0]
    methane = AdvancedThermo.get_fluid_component("methane")
    m = AdvancedThermo.create_mixture([methane])
    da = AdvancedThermo.derivatives_of_a(m,AdvancedThermo.ρ(m,p,T,x),T,x)
    
    pc = 5.0e6
    Tc = 192.0
    ρc = 10.0e3
    methane_new = AdvancedThermo.get_fluid_component("methane")
    methane_new = AdvancedThermo.new_critical!(methane_new,pc,Tc,ρc)
    m_new = AdvancedThermo.create_mixture([methane_new])
    da_new = AdvancedThermo.derivatives_of_a(m_new,AdvancedThermo.ρ(m_new,p,T,x),T,x)
    
    # compare Helmholtz energy and derivatives (new vs. old)
    tol = 5.e-13
    val = DiffResults.value(da)
    val_new = DiffResults.value(da_new)
    grad = DiffResults.gradient(da)
    grad_new = DiffResults.gradient(da_new)
    hess = DiffResults.hessian(da)
    hess_new = DiffResults.hessian(da_new)
    @test abs((val-val_new)/val)<tol
    @test abs((grad[1]-grad_new[1])/grad[1])<tol
    @test abs((grad[2]-grad_new[2])/grad[2])<tol
    @test abs((grad[3]-grad_new[3])/grad[3])<tol
    @test abs((hess[1,1]-hess_new[1,1])/hess[1,1])<tol
    @test abs((hess[1,2]-hess_new[1,2])/hess[1,2])<tol
    @test abs((hess[1,3]-hess_new[1,3])/hess[1,3])<tol
    @test abs((hess[2,2]-hess_new[2,2])/hess[2,2])<tol
    @test abs((hess[2,3]-hess_new[2,3])/hess[2,3])<tol
    @test abs((hess[3,3]-hess_new[3,3])/hess[3,3])<tol
end

I call this a great success! The code so far can be obtained by downloading the zip file. Now that I have a proven way to modify the critical parameters for a pure fluid, the next part of the series will examine the effect on mixtures with helium.

References

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

[2] Messerly et al., J. Chem. Eng. Data, 2020, 65, 1028 (link to publisher)

[3] 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