Advanced Thermodynamics

Thermodynamic properties of natural gas mixtures containing helium - Part 9

7 August 2021

This is Part 9 of a series.

In Part 8, I validated my earlier calculations for mixtures with helium against TREND 4.0 [1] and trialled a different mixing rule for the binary interaction parameters between methane and helium, with improved results for density predictions. In this part, I'll see whether the predictions of phase equilibria also improve.

Using the arithmetic_γ function developed in Part 8, it's a simple matter to compare against the pTxy dataset as was done in Part 7 using the 'classical' critical constants for helium reported by Messerly et al. [2]. The important statistics are given below:

    Results for γv = γT = 1
    -----------------------
    Statistic      lnK(methane)               lnK(helium)
    AAD            0.39529678667060447        8.089845336234148
    Bias           0.39393776759098836        -8.089845336234148
    RMS            0.5033200776369687         15.277498663484204
    Max            1.9302763783863863         -109.38993501770257

    Results for γv = 1.086694, γT = 2.044932
    ----------------------------------------
    Statistic      lnK(methane)               lnK(helium)
    AAD            0.2141818363559133         4.401756605836159
    Bias           0.2092722154474173         -4.401756605836159
    RMS            0.3363407531427379         9.216561743227535
    Max            1.8747999120906926         -67.62638550447191

I haven't yet looked in detail at the cause for the large Max deviations. With the GERG-2008 EOS [3], a point where my flash leads to anomalously large lnK values is:

    "Reference",            "n", "T",     "p",     "x1",    "x2",    "y1",    "y2"
    "DeVaney et al., 1971", 10,  154.01,  2757.9,  0.9950,  0.0050,  0.5258,  0.4742

One possibility for creating a more robust solver is to use a sophisticated technique like in IntervalRootFinding.jl, the documentation for which can be found here. The phase equilibrium problem is fairly technical to try first, so instead I'll try making a more robust density solver than the one in Part 4. Setting up the call to IntervalRootFinding is straightforward:

function ρ_robust(m::FluidMixture,p_target,T,x)
    return IntervalRootFinding.roots(ρ->p(m,ρ,T,x)-p_target,1.e-4..1.e5)
end

@time AdvancedThermo.ρ_robust(m,1.e6,T,x)
    2.209085 seconds (6.28 M allocations: 244.578 MiB, 5.04% gc time)
  1-element Vector{Root{Interval{Float64}}}:
   Root([736.006, 736.007], :unique)

@time AdvancedThermo.ρ_robust(m,3.3e6,T,x)
   23.685067 seconds (68.19 M allocations: 2.595 GiB, 5.10% gc time)
  3-element Vector{Root{Interval{Float64}}}:
   Root([3873.39, 3873.4], :unique)
   Root([17234.3, 17234.4], :unique)
   Root([7522.07, 7522.08], :unique)

Last time, I set up TREND [1] in a way that allowed me to calculate the properties of mixtures containing helium using the modified critical parameters from Messerly et al. [2]. As TREND has a phase equilibrium solver, I can use its Excel API to calculate the flash results. What I find is somewhat encouraging but not conclusive:

• When I flash the points in the file pTxy_methane_helium_pruned.csv with the GERG-2008 EOS [3], TREND 4.0 reports a single-phase result for almost 20% of cases, and at least one point leads to an infinite loop
• With the model using Messerly et al. critical parameters and γ = 1, all points are reported as two-phase but the predicted K_helium are usually several orders of magnitude larger than experiment
• With the model using Messerly et al. critical parameters, γ_v = 1.088694 and γ_T = 2.044932, all points are two-phase. The predicted K_helium are still several orders of magnitude larger than experiment (but better than γ = 1)

As TREND failed to find the correct two-phase equilibrium in a large percentage of cases with the GERG-2008 EOS, I can't compute a full complement of statistics. So I really need to find a more robust way to flash these mixtures before I can say which model is better. I happen to also have REFPROP 10 [4] available: it uses the same fluid files as TREND so, in principle, it should be relatively straightforward to try out the new model with. Unfortunately, when I try to flash the mixtures in pTxy_methane_helium.csv with the GERG-2008 EOS from REFPROP 10:

• 21 points are found to be single phase
• another point has phase fraction around 4×10⁻¹²
• 1 point fails in the density solver
• 2 points fail to converge in the flash

Discussion

At this point, I've exhausted the external tools available to me. Although they were not completely up-to-task, I was able to make some progress. My sense is that some set of critical parameters for helium will improve the predictive capability of the EOS, but that it may not be those of Messerly et al. [2], i.e. the critical temperature and critical density may need to be optimised in a similar way as we did when using the cubic EOS [5]. Either way, to reach a definitive conclusion, it will be necessary to have a more robust phase equilibrium solver.

Building a robust phase equilibrium solver is not an easy or straightforward task; I would have preferred to find an accessible tool to perform the calculations for me in a reliable way. It's disappointing to have investigated this far only to come up against a sizeable obstacle. However, I believe it is surmountable. I will pursue the development of a robust solver in another (hopefully short) series, and will resume this series once I have the necessary tools.

The code so far can be obtained by downloading the zip file.

References

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

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

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

[4] Lemmon et al., REFPROP Documentation, Release 10.0, NIST (link to FAQ)

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