Advanced Thermodynamics

Thermodynamic properties of natural gas mixtures containing helium - Part 1

19 June 2021

This piece is Part 1 of a series.

In 2017, Tom Hughes, Eric May and I published the paper "Effective Critical Constants for Helium for Use in Equations of State for Natural Gas Mixtures" [1]. In that paper we used the Peng-Robinson cubic equation of state and modified the critical temperature and pressure of helium to improve predictions of vapour-liquid equilibrium data (compared to the case where the true critical constants are used). Near the end of that paper we wrote

"Improving property correlations and predictions by multiparameter Helmholtz equations of state for helium-containing fluids is an avenue for future work."

To begin, I need an implementation of a multiparameter Helmholtz equation of state (EOS). As I potentially want to make modifications to various parts of the EOS I'll code the necessary parts of the EOS from scratch, based on the GERG-2008 EOS developed by Kunz and Wagner [2]. The GERG-2008 EOS is nice-to-use for this purpose as the reduced residual Helmholtz energy α^r of each pure fluid is described by a straightforward summation of 'polynomial' and 'exponential' functions of the reduced density δ and inverse reduced temperature τ, viz.

The code will be written in Julia. Hopefully the code is simple enough to understand even if this is the first time you are seeing Julia but you can find many guides online if required.

In Julia, the functions for summing the polynomial terms and exponential terms can be written succinctly like the following:

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

mutable struct FluidComponent
  Name::String
  Tc::Float64   # critical temperature [K]
  Pc::Float64   # critical pressure [Pa]
  ρc::Float64   # critical density [mol/m^3]
  mw::Float64   # molar mass [kg/mol]
  # polynomial-type parameters
  np::Vector{Float64}
  dp::Vector{Int64}
  tp::Vector{Float64}
  # exponential-type parameters
  ne::Vector{Float64}
  de::Vector{Int64}
  te::Vector{Float64}
  ce::Vector{Int64}
end

function αʳ(δ,τ,f::FluidComponent)
  return αʳ(δ,τ,f.np,f.dp,f.tp) + αʳ(δ,τ,f.ne,f.de,f.te,f.ce)
end

At this point I need a fluid to work with. As helium will shortly be the focus, I'll create a function that prepares a FluidComponent for helium [2]:

function get_fluid_component(name)
  if name=="helium"
      # refer Kunz and Wagner (2012) Tables A4-A5
      np = [
          -0.45579024006737
          0.12516390754925e1
          -0.15438231650621e1
          0.20467489707221e-1
      ]
      dp = [
          1
          1
          1
          4
      ]
      tp = [
          0.000
          0.125
          0.750
          1.000
      ]
      ne = [
          -0.34476212380781
          -0.20858459512787e-1
          0.16227414711778e-1
          -0.57471818200892e-1
          0.19462416430715e-1
          -0.33295680123020e-1
          -0.10863577372367e-1
          -0.22173365245954e-1
      ]
      de = [
          1
          3
          5
          5
          5
          2
          1
          2
      ]
      te = [
          0.750
          2.625
          0.125
          1.250
          2.000
          1.000
          4.500
          5.000
      ]
      ce = [
          1
          1
          1
          1
          1
          2
          3
          3
      ]
      return FluidComponent(
          name,
          5.1953,
          0.0, # NOTE(DR) critical pressure not required (yet)
          17399.0,
          4.002602e-3,
          np,
          dp,
          tp,
          ne,
          de,
          te,
          ce
      )
  else
      # TODO(DR): more fluids as needed
  end
end

function αʳ(fluid::FluidComponent,ρ,T)
  δ = ρ/fluid.ρc
  τ = fluid.Tc/T
  return αʳ(δ,τ,fluid)
end

To test the code, I'll use TREND 4.0 [3] to make some calculations of the reduced Helmholtz energy of pure helium with the GERG-2008 EOS [2] and compare the results.

Table 1. Comparison between the current EOS implementation for helium and results from the GERG-2008 implementation in TREND 4.0.

@testset "αʳ(helium)" begin
  helium = AdvancedThermo.get_fluid_component("helium")
  reltol = 1.e-14
  # benchmark values from TREND 4.0 (GERG-2008)
  αʳ = AdvancedThermo.αʳ(helium,50.0,300.0)
  @test abs(αʳ-5.97611193523994E-04)/αʳ<reltol
  αʳ = AdvancedThermo.αʳ(helium,20000.0,100.0)
  @test abs(αʳ-2.709103301953120E-01)/αʳ<reltol
end

  Testing Running tests...
  Test Summary: | Pass  Total
  αʳ(helium)    |    2      2
       Testing AdvancedThermo tests passed

The code so far can be obtained by downloading the zip file. Upcoming articles will build on this work and (hopefully) lead to some novel results. Click here for Part 2.

References

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

[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