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."This post is a first step in that direction.

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.*

Polynomial terms: | ∑ nδ^{d}τ^{t} |

Exponential terms: | ∑ nδexp(−^{d}τ^{t}δ)^{c} |

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
```

To make it easier to keep the parameters for a particular fluid together I'll
create a struct to contain them like so:
```
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
```

then the reduced Helmholtz energy of a `FluidComponent`

will be calculated with the following function:
```
function αʳ(δ,τ,f::FluidComponent)
return αʳ(δ,τ,f.np,f.dp,f.tp) + αʳ(δ,τ,f.ne,f.de,f.te,f.ce)
end
```

Note that I've used the same name (`αʳ`

) for several functions: this
is fine in Julia as long as the argument lists are distinctive.
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
```

The final piece of code I'd like at this stage is a function
that accepts the absolute density and temperature, instead of their
reduced versions. This can be written like so:
```
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.

ρ /mol·m^{−3} | T /K |
α^{r} (this work) | α^{r} [3] |
---|---|---|---|

50.0 | 300.0 | 5.976111935239946×10^{−4} |
5.97611193523994×10^{−4} |

20000.0 | 100.0 | 0.2709103301953123 | 0.270910330195312 |

```
@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
```

Running the tests gives the following output:
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.

[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

Go to Home | How to cite this page | Permalink to this page |