Advanced Thermodynamics

Thermodynamic properties of natural gas mixtures containing helium - Part 3

3 July 2021

This is Part 3 of a series.

In Part 2, I showed how to calculate the residual contribution to the reduced Helmholtz energy of mixtures based on the mixing rules developed by Kunz et al. [1] as part of the GERG-2004 EOS. To calculate other fluid properties — such as pressure, enthalpy, or speed of sound — derivatives of the Helmholtz energy are needed. The task of calculating derivatives will be explored below.

One of the main reasons that I chose Julia to implement the Equation of State was that some Automatic Differentiation (AD) packages are available. Although the derivatives could be coded by hand, this will hopefully provide a shortcut. There are actually a large number of available packages: I'll try ForwardDiff.jl at first.

Relationships between the derivatives of the Helmholtz energy and common thermodynamic properties are given in Table B1 of Kunz and Wagner [2]. There are two main things it looks like I need to address:

1. the derivatives are in terms of Helmholtz energy a, not α^r
2. derivatives are with respect to v, T and x

function αᵒ(fluid::FluidComponent,ρ,T)
    δᵢ = ρ/fluid.ρc
    τᵢ = fluid.Tc/T
    s = [
        1.0, 
        τᵢ, 
        log(τᵢ), 
        log(abs(sinh(fluid.vi[4]*τᵢ))), 
       -log(abs(cosh(fluid.vi[5]*τᵢ))), 
        log(abs(sinh(fluid.vi[6]*τᵢ))), 
       -log(abs(cosh(fluid.vi[7]*τᵢ)))]
    return log(δᵢ) + fluid.ni's
end
function αᵒ(m::FluidMixture,ρ,T,x)
    result = 0.0
    for i = 1:length(x)
        if x[i]>0.0
            result += x[i]*(αᵒ(m.components[i],ρ,T)+log(x[i]))
        end
    end
    return result
end

    # ideal-gas terms
    ni::Vector{Float64}
    vi::Vector{Float64}

The reduced Helmholtz energy is the sum of the ideal-gas contribution and the residual contribution:

function α(m::FluidMixture,ρ,T,x)
    return αᵒ(m,ρ,T,x) + αʳ(m,ρ,T,x)
end

Regarding point 2, from Table B1 [2] the properties that will be most readily calculable are the following:

• p = −(∂a/∂v)_T,x
• s = −(∂a/∂T)_v,x

• u = a + Ts
• h = u + pv
• g = a + pv

function a(m::FluidMixture,v,T,x)
  return R*T*α(m,inv(v),T,x)
end
function ∂a∂T(m::FluidMixture,v,T,x)
  f = T -> a(m,v,T,x)
  return ForwardDiff.derivative(f,T)
end
function s(m::FluidMixture,ρ,T,x)
  return -∂a∂T(m,inv(ρ),T,x)
end

Once again, I'll compare with TREND 4.0 [3] and create a unit test:

@testset "entropy (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)
  s = AdvancedThermo.s(mixture,1000.0,300.0,x)
  @test abs(s-(-21.682667884))<1.e-3
end

function ∂a∂v(m::FluidMixture,v,T,x)
  f = v -> a(m,v,T,x)
  return ForwardDiff.derivative(f,v)
end
function p(m::FluidMixture,ρ,T,x)
  return -∂a∂v(m,inv(ρ),T,x)
end

@testset "pressure (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)
  p = AdvancedThermo.p(mixture,1000.0,300.0,x)
  @test abs(p-2.514507019437e6)<5.e-3
end

I could continue bulding up the API in the same manner as above to create functions for g, u, h etc. However, when reading over the ForwardDiff advanced documentation it looks like it will be beneficial to use DiffResults.jl in conjunction with ForwardDiff.jl. The documentation at the time I'm writing this can be found here. For example, code to evaluate all derivatives of the Helmholtz energy up to second order with respect to volume, temperature and mole fraction composition looks like:

function derivatives_of_a(m::FluidMixture,ρ,T,x)
  z = vcat(inv(ρ),T,x)
  result = DiffResults.HessianResult(z)
  f = z -> a(m,z[1],z[2],z[3:end])
  result = ForwardDiff.hessian!(result, f, z)
  return result
end

The result from DiffResults.jl is returned in an object of type DiffResults.MutableDiffResult. A better API than what I have implemented so far will be needed. However for now, I have added a unit test that verifies some of the first- and second-order derivatives against TREND 4.0 [3] and serves to demonstrate the functionality:

@testset "derivatives of a" begin
  nitrogen = AdvancedThermo.get_fluid_component("nitrogen")
  helium = AdvancedThermo.get_fluid_component("helium")
  mixture = AdvancedThermo.create_mixture([nitrogen,helium])
  x = [0.7, 0.3]
  result = AdvancedThermo.derivatives_of_a(mixture,1000.0,300.0,x)
  # Benchmark values from TREND 4.0
  # helmholtz energy a
  @test abs(DiffResults.value(result)-3991.1843584)<0.1
  # p = -∂a/∂v
  @test abs(-DiffResults.gradient(result)[1]-2.514507019437e6)<5.e-3
  # s = -∂a/∂T
  @test abs(-DiffResults.gradient(result)[2]-(-21.682667884))<5.e-3
  # ∂p/∂T_v = -∂²a/∂v∂T
  @test abs(-DiffResults.hessian(result)[1,2]-0.00860851125970e6)<1.e-4
  # volume expansivity = 1/v * ∂v/∂T_p = 1/v * -(∂p/∂T)_v/(∂p/∂v)_T
  v_exp = -1000.0*DiffResults.hessian(result)[1,2]/DiffResults.hessian(result)[1,1]
  @test abs(v_exp-0.003392906770418)<1.e-8 
end

Of the quantities that aren't simplistically related to derivatives of a, the partial fugacity coefficients are perhaps the most important for examining phase equilibria, which will be one of the next phases of this series. The relation in [1] for the partial fugacity coefficient of component i is

• lnφ_i = (∂nα^r/∂n_i) _T,V,n' − lnZ

function nαʳ(m::FluidMixture,VTn)
  ntotal = sum(VTn[3:end])
  return ntotal*αʳ(m,ntotal/VTn[1],VTn[2],VTn[3:end]./ntotal)
end
function derivatives_of_nαʳ(m::FluidMixture,V,T,n)
    z = vcat(V,T,n)
    result = DiffResults.HessianResult(z)
    f = z -> nαʳ(m,z)
    result = ForwardDiff.hessian!(result, f, z)
    return result
end

@testset "derivatives of nαʳ" begin
  nitrogen = AdvancedThermo.get_fluid_component("nitrogen")
  helium = AdvancedThermo.get_fluid_component("helium")
  mixture = AdvancedThermo.create_mixture([nitrogen,helium])
  x = [0.7, 0.3]
  ρ = 1000.0
  T = 300.0
  result = AdvancedThermo.derivatives_of_nαʳ(mixture,inv(ρ),T,x)
  Z = AdvancedThermo.Z(mixture,ρ,T,x)
  ϕ_N2 = exp(DiffResults.gradient(result)[3] - log(Z))
  ϕ_He = exp(DiffResults.gradient(result)[4] - log(Z))
  @test abs(ϕ_N2-0.998810946)<1.e-9
  @test abs(ϕ_He-1.028530046)<1.e-9
end

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] Kunz and Wagner, J. Chem. Eng. Data, 2012, 57, 3032 (link to publisher)

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