Advanced Thermodynamics

Thermodynamic properties of natural gas mixtures containing helium - Part 4

4 July 2021

This is Part 4 of a series.

In Part 3, I showed how to calculate derivatives of the Helmholtz energy in Julia via Automatic Differentiation. In this part, a density solver will be implemented so that I can specify pressure and temperature instead of density and temperature, which will help when examining experimental data.

A consequence of using multi-parameter equations of state in the GERG-2008 model is that multiple density solutions/roots can exist at a particular pressure and temperature [1]. However, as noted by Kunz et al.:

"the multiple density solutions can be handled in a quite efficient and reliable way, without any problems concerning the calculation of thermodynamic properties for the various types of mixtures investigated in this work"

• for the vapour region, start at low density and search upward in density until either
  • the specified pressure is found
  • a local maximum is found (therefore the vapour root does not exist)
• for the liquid region, start at high density and search downward until either
  • the specified pressure is found
  • a local minimum is found (therefore the liquid root does not exist)
• if two roots are found, return the one with lower Gibbs energy

function ρ(m::FluidMixture,p,T,x)
  ρ_lo = 1.e-6
  (p_lo,∂p∂ρ,∂²p∂ρ²) = derivatives_of_p(m,ρ_lo,T,x)
  while abs(p_lo-p)>1.e-4
      dρ = (p-p_lo)/∂p∂ρ
      if ρ_lo + dρ > 0.0
          ρ_lo += dρ
          (p_lo,∂p∂ρ,∂²p∂ρ²) = derivatives_of_p(m,ρ_lo,T,x)
      else
          ρ_lo = -1.0
          p_lo = p
      end
  end
  ρ_hi = 1.e6
  (p_hi,∂p∂ρ,∂²p∂ρ²) = derivatives_of_p(m,ρ_hi,T,x)
  while abs(p_hi-p)>1.e-4
      dρ = (p-p_hi)/∂p∂ρ
      dρ_s = -∂p∂ρ/∂²p∂ρ²
      if dρ>dρ_s
          ρ_hi += dρ
          (p_hi,∂p∂ρ,∂²p∂ρ²) = derivatives_of_p(m,ρ_hi,T,x)
      else
          ρ_hi = -1.0
          p_hi = p
      end
  end
  if ρ_lo>=0.0 && ρ_hi>=0.0
      g_lo = g(m,ρ_lo,T,x)
      g_hi = g(m,ρ_hi,T,x)
      if g_lo < g_hi
          return ρ_lo
      else
          return ρ_hi
      end
  elseif ρ_lo >= 0.0
      return ρ_lo
  else
      return ρ_hi
  end
end

function g(m::FluidMixture,ρ,T,x)
  v = inv(ρ)
  return a(m,v,T,x) + p(m,ρ,T,x)*v
end

function derivatives_of_p(m::FluidMixture,ρ,T,x)
  result = DiffResults.HessianResult([ρ])
  f = ρ -> p(m,ρ[1],T,x)
  result = ForwardDiff.hessian!(result, f, [ρ])
  pcalc = DiffResults.value(result)
  ∂p∂ρ = DiffResults.gradient(result)[1]
  ∂²p∂ρ² = DiffResults.hessian(result)[1,1]
  return (pcalc,∂p∂ρ,∂²p∂ρ²)
end

@testset "gibbs energy (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)
  g = AdvancedThermo.g(mixture,1000.0,300.0,x)
  @test abs(g-6505.6913733708)<1.e-1
end
@testset "density solver (nitrogen)" begin
  nitrogen = AdvancedThermo.get_fluid_component("nitrogen")
  mixture = AdvancedThermo.create_mixture([nitrogen])
  # benchmark values from TREND 4.0 (GERG-2008)
  @test abs(AdvancedThermo.ρ(mixture,1.e4,100.0,[1.0])-12.0505465245949)<1.e-7
  @test abs(AdvancedThermo.ρ(mixture,1.e6,100.0,[1.0])-24661.1684674037)<1.e-7
  @test abs(AdvancedThermo.ρ(mixture,1.e6,115.0,[1.0])-1223.49168349671)<1.e-7
  @test abs(AdvancedThermo.ρ(mixture,1.e7,115.0,[1.0])-23571.4077503479)<1.e-7
  @test abs(AdvancedThermo.ρ(mixture,2.e6,126.0,[1.0])-2500.28432820219)<1.e-7
  @test abs(AdvancedThermo.ρ(mixture,3.39e6,126.0,[1.0])-14332.0445222086)<2.e-6
  @test abs(AdvancedThermo.ρ(mixture,1.e7,126.0,[1.0])-21321.8153441504)<1.e-7

end
@testset "density solver (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)
  @test abs(AdvancedThermo.ρ(mixture,1.e4,300.0,x)-4.0089604347598)<1.e-6
  @test abs(AdvancedThermo.ρ(mixture,1.e4,100.0,x)-12.0356259880732)<1.e-8
  @test abs(AdvancedThermo.ρ(mixture,1.e9,300.0,x)-46376.8316627108)<1.e-8
end

I want to work on the phase equilibrium solver in the next part. In anticipation of that, I've implemented a simple function for calculating the partial fugacity coefficients (as natural log) of all components in the FluidMixture:

function lnϕ(m::FluidMixture,ρ,T,x)
  result = DiffResults.GradientResult(x)
  f = x -> nαʳ(m,vcat(inv(ρ),T,x))
  result = ForwardDiff.gradient!(result, f, x)
  return DiffResults.gradient(result) .- log(Z(m,ρ,T,x))
end

  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

  lnϕ = AdvancedThermo.lnϕ(mixture,ρ,T,x)
  @test abs(exp(lnϕ[1])-0.998810946)<1.e-9
  @test abs(exp(lnϕ[2])-1.028530046)<1.e-9

The code so far can be obtained by downloading the zip file. In the next part of the series, I'll implement a simple phase equilibrium solver: soon it will be possible to examine experimental pTxy data.

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] Topliss, Techniques to Facilitate the Use of Equations of State for Complex Fluid-Phase Equilibria, PhD Dissertation, University of California, Berkeley, 1985

[4] Bell and Alpert, Fluid Phase Equilib., 2018, 476, 89 (link to publisher)