Additional properties of a mixture

[infinity77]

Hi,

I have tried using phasepy to analyze the behavior of a mixture of 11 components - as a gas - and comparing the results with the AGA8 code (https://pages.nist.gov/AGA8/). I have been looking specifically at density, Z, entropy, enthalpy, heat capacities and speed of sound.

Quite curiously, I get very similar results for density and Z between phasepy and AGA8, but all the other properties are very far off from one another. This is what I get:

Property	PhasePy	AGA8
Molar density (mol/l)	0.0424631	0.0424358
Compressibility factor (Z)	0.996037	0.996618
Residual Entropy (J/mol-K)	-0.0671334	5.050331
Residual Enthalpy (J/mol)	-28.8389	-439.650
Residual isochoric heat capacity (J/mol-K)	0.017577	32.6788
Residual isobaric heat capacity (J/mol-K)	0.152745	41.1268
Speed of sound (m/s)	437.46	379.25

I have verified the AGA8 results with the Excel interface and the Fortran interface. The code I have used for phasepy is below:

from __future__ import division, print_function, absolute_import

import numpy as np
from phasepy import preos, component, mixture

# KIJ interaction coefficients
#       Methane,  Nitrogen  CO2     Ethane   Propane Isobutane Butane Isopentane Pentane  Hexane   H2S
KIJ = [[0.00000, 0.02890, 0.09780, -0.0059, 0.01190, 0.02560, 0.01850, -0.0056, 0.02300, 0.04000, 0.00000],
       [0.02890, 0.00000, -0.0122, 0.05330, 0.08780, 0.10330, 0.07110, 0.09220, 0.10000, 0.14960, 0.16520],
       [0.09780, -0.0122, 0.00000, 0.13000, 0.13150, 0.13000, 0.13520, 0.12190, 0.12520, 0.11000, 0.09670],
       [-0.0059, 0.05330, 0.13000, 0.00000, 0.00110, -0.0067, 0.00890, 0.00000, 0.00780, -0.0400, 0.09520],
       [0.01190, 0.08780, 0.13150, 0.00110, 0.00000, -0.0078, 0.00330, 0.01110, 0.02670, 0.00070, 0.08780],
       [0.02560, 0.10330, 0.13000, -0.0067, -0.0078, 0.00000, -0.0004, 0.00000, 0.00000, 0.00000, 0.04740],
       [0.01850, 0.07110, 0.13520, 0.00890, 0.00330, -0.0004, 0.00000, 0.00000, 0.01740, -0.0056, 0.00000],
       [-0.0056, 0.09220, 0.12190, 0.00000, 0.01110, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000],
       [0.02300, 0.10000, 0.12520, 0.00780, 0.02670, 0.00000, 0.01740, 0.00000, 0.00000, 0.00000, 0.06300],
       [0.04000, 0.14960, 0.11000, -0.0400, 0.00070, 0.00000, -0.0056, 0.00000, 0.00000, 0.00000, 0.00000],
       [0.00000, 0.16520, 0.09670, 0.09520, 0.08780, 0.04740, 0.00000, 0.00000, 0.06300, 0.00000, 0.00000]]

KIJ = np.array(KIJ)

# Critical properties of pure chemicals
#        Mw       Tc (K)  Pc (bar) Vc(cm3/mol) Zc     omega
DATA = [[16.043 , 190.56 , 46.00  , 99.0     , 0.288, 0.011],
        [28.0135, 126.2  , 33.94  , 89.5     , 0.29 , 0.040],
        [44.01  , 304.2  , 73.76  , 94.0     , 0.274, 0.228],
        [30.07  , 305.36 , 48.80  , 146.0    , 0.284, 0.099],
        [44.097 , 369.9  , 42.50  , 199.0    , 0.281, 0.152],
        [58.123 , 407.84 , 36.40  , 256.0    , 0.282, 0.177],
        [58.123 , 425.2  , 37.90  , 257.0    , 0.274, 0.193],
        [72.15  , 460.37 , 33.50  , 313.0    , 0.27 , 0.228],
        [72.15  , 469.7  , 33.70  , 310.0    , 0.269, 0.249],
        [86.177 , 507.5  , 30.30  , 366.0    , 0.264, 0.305],
        [34.082 , 373.2  , 89.37  , 98.5     , 0.284, 0.245]]

DATA = np.array(DATA)

data = [('Methane', 0.8256), ('Nitrogen', 0.0116) , ('CO2', 0.0135)   , ('Ethane', 0.0572), 
        ('Propane', 0.0449), ('Isobutane', 0.0117), ('Butane', 0.0195), ('Isopentane', 0.0085), 
        ('Pentane', 0.0038), ('Hexane', 0.0025)   , ('H2S', 0.0012)]

row = 0
phasepy_components = []
compositions = []

for name, fraction in data:
    comp = component(name=name, Tc=DATA[row, 1], Pc=DATA[row, 2], Vc=DATA[row, 3],
                     Zc=DATA[row, 4], w=DATA[row, 5], Mw=DATA[row, 0])

    phasepy_components.append(comp)
    compositions.append(fraction)

    if row == 1:
        # Create th mixture as 2 components are done already
        mix = mixture(phasepy_components[0], phasepy_components[1])
    elif row > 1:
        # Add to the existing mixture
        mix.add_component(phasepy_components[-1])

    row += 1

mix.kij_cubic(KIJ)
eos = preos(mix, mixrule='qmr')

T = 288.15   # K
P = 1.01325  # bar
x = np.array(compositions)
print('Molar density                    : ', 1e3*eos.density(x, T, P, 'V'), 'mol/l')
print('Compressibility factor (Z)       : ', eos.Zmix(x, T, P)[0])

# Thermal derived properties
Sr = eos.EntropyR(x, T, P, 'V')
Hr = eos.EnthalpyR(x, T, P, 'V')
Cvr = eos.CvR(x, T, P, 'V')
Cpr = eos.CpR(x, T, P, 'V')

# ideal gas heat capacities, better values can be obtained with DIPPR 801 correlations
r = 8.314  # J / mol K
CvId = 3*r/2
CpId = 5*r/2
w = eos.speed_sound(x, T, P, 'V', CvId=CvId, CpId=CpId)

print('Residual Entropy                 : ', Sr , 'J/mol-K')
print('Residual Enthalpy                : ', Hr , 'J/mol')
print('Residual isochoric heat capacity : ', Cvr, 'J/mol-K')
print('Residual isobaric heat capacity  : ', Cpr, 'J/mol-K')
print('Speed of sound                   : ', w  , 'm/s')

I am sure I am missing something obvious - could you please explain what that would be?

Thank you in advance.

[tkeskita]

Hi, can you share your AGA8 Fortran code? Are you sure its giving you the residual properties?

[gustavochm]

Hi!

Thanks for you interest in Phasepy.

The thermal derivatives properties (Entropy, Enthalpy, Cv, Cp and speed of sound) were just added recently to the package. After this enquiry I checked the functions and there is no bug in the functions (the relations are correctly coded). I used the following expressions:

There are several reasons that could explain the differences between the AGA8 package.

• In phasepy the thermal derivatives of the Helmholtz free energy are computed numerically (I tried different steps at your code and this didn't changed the results)
• Could you please confirm that the reported values from the Entropy and Enthalpy from AGA8 are computed in the (T,P,N) ensemble and that the reported value is the residual one.
• For the heat capacities it is my belief that the AGA8 package might be returning the total heat capacities (residual + ideal). Consider that the ideal contribution takes the values reported for methane (CpId = 35.315 # J / mol K and CvId = CpId - r # J / mol K, obtained from DIPPR 801), then the computed total isochoric and isobaric heat capacities are 27.01 and 35.46 J/mol K , respectively. Way closer to the reported values by AGA8. To complement this point, the reported isochoric and isobaric heat capacities on NIST for methane at 288K and 1 bar are 27.050 and 35.448 J / mol K.
• For the speed of sound calculation, the result can change considerably by means of the ideal gas contribution (CvId, CpId). Better results of this contribution can be obtained from the DIPPR 801 correlation. For example, as your mixture is mostly methane its reported value could be used (CpId = 35.315 # J / mol K and CvId = CpId - r # J / mol K). With those values the computed speed of sound is 387.097 m/s. This result is closer to the reported by AGA8.

It would be interesting to compare the source codes of both implementations!

I hope this help!

[infinity77]

Hi,

you are indeed correct and the AGA8 code is not reporting residual properties, I completely missed it from looking at the code.

the AGA8 implementation (at least for simple stuff) is public domain and you can find it here in 4 different languages:

https://github.com/usnistgov/AGA8