Issue with TTSE interpolation when using setState_pT for supercritical CO2 about externalmedia

[casella]

It seems that setState_pT picks the wrong enthalpy for some reason:

[casella]

Unfortunately this issue is still there in version 3.3.2, updated to the latest version of CoolProp.

[bilderbuchi]

This is maybe an upstream (i.e. Coolprop) issue, I can reproduce this reading using the Python wrapper:

In [1]: import CoolProp

In [2]: HEOS = CoolProp.AbstractState("HEOS", "CO2")

In [3]: TTSE = CoolProp.AbstractState("TTSE&HEOS", "CO2")

In [4]: BICU = CoolProp.AbstractState("BICUBIC&HEOS", "CO2")

In [5]: HEOS.update(CoolProp.PT_INPUTS, 90e5, 308); BICU.update(CoolProp.PT_INPUTS, 90e5, 308); TTSE.update(CoolProp.PT
   ...: _INPUTS, 90e5, 308)

In [7]: print(HEOS.rhomass(), TTSE.rhomass(), BICU.rhomass())
665.3674051097354 905.0585107989239 636.4883422212183

In [9]: CoolProp.__version__
Out[9]: '6.4.3'

In [10]: HEOS.update(CoolProp.PT_INPUTS, 90e5, 306); BICU.update(CoolProp.PT_INPUTS, 90e5, 306); TTSE.update(CoolProp.P
    ...: T_INPUTS, 90e5, 306)

In [11]: print(HEOS.rhomass(), TTSE.rhomass(), BICU.rhomass())
702.8534640782832 705.8648349168404 684.6191073389101

(checked two relevant points only, irrelevant lines elided)

I looked in the issue tracker, the most relevant-sounding issues I could find are CoolProp/CoolProp#1301 CoolProp/CoolProp#1437 (for the bicubic method, though)

@jowr ?

[bilderbuchi]

I have adapted the script from here to show the issue. This seems to show some periodicity, but especially severe in your region of interest:

Graph generation with Python wrapper

import CoolProp
import CoolProp.CoolProp as CP
import matplotlib.pyplot as plt
import matplotlib.colors as colors
import matplotlib.cm as cmx
import matplotlib.ticker
import numpy as np
import random

fig = plt.figure(figsize=(10, 5))
ax1 = fig.add_axes((0.08, 0.1, 0.32, 0.83))
ax2 = fig.add_axes((0.50, 0.1, 0.32, 0.83))

Ref = 'CO2'
p_lim = (70e5, 100e5)
T_lim = (290, 330)
# adjust limits of error color bar in definition of cNorm
cmap = 'plasma'

BICUBIC = CoolProp.AbstractState('BICUBIC&HEOS', Ref)
TTSE = CoolProp.AbstractState('TTSE&HEOS', Ref)
EOS = CoolProp.AbstractState('HEOS', Ref)

Ts = []
ps = []
for T_trial in np.linspace(T_lim[0], T_lim[1], 1000):
    try:
        p_sat = CP.PropsSI('P', 'T', T_trial, 'Q', 0, Ref)
        ps.append(p_sat)
        Ts.append(T_trial)
    except ValueError as exc:
        pass  # Temperature beyond Tcrit

TTT1, HHH1, PPP1, EEE1 = [], [], [], []
TTT2, HHH2, PPP2, EEE2 = [], [], [], []

cNorm = colors.LogNorm(vmin=1e-3, vmax=100)
scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=plt.get_cmap(cmap))

for a_useless_counter in range(40000):

    # h = random.uniform(*h_lim)
    T = random.uniform(*T_lim)
    p = 10**random.uniform(np.log10(p_lim[0]), np.log10(p_lim[1]))
    CP.set_debug_level(0)
    try:

        EOS.update(CoolProp.PT_INPUTS, p, T)
        rhoEOS = EOS.rhomolar()
        hEOS = EOS.hmass()

        TTSE.update(CoolProp.PT_INPUTS, p, T)
        rhoTTSE = TTSE.rhomolar()
        hTTSE = TTSE.hmass()

        BICUBIC.update(CoolProp.PT_INPUTS, p, T)
        rhoBICUBIC = BICUBIC.rhomolar()
        hBICUBIC = BICUBIC.hmass()

        errorTTSE = abs(rhoTTSE/rhoEOS-1)*100
        errorBICUBIC = abs(rhoBICUBIC/rhoEOS-1)*100
        if errorTTSE > 100 or errorTTSE < 1e-12:
            print(T, p, errorTTSE)

        TTT1.append(T)
        HHH1.append(hTTSE)
        PPP1.append(p)
        EEE1.append(errorTTSE)

        TTT2.append(T)
        HHH2.append(hBICUBIC)
        PPP2.append(p)
        EEE2.append(errorBICUBIC)

    except ValueError as VE:
        print('ERROR', VE)
        pass

SC1 = ax1.scatter(TTT1, PPP1, s=8, c=EEE1, edgecolors='none',
                  cmap=plt.get_cmap(cmap), norm=cNorm)
SC2 = ax2.scatter(TTT2, PPP2, s=8, c=EEE2, edgecolors='none',
                  cmap=plt.get_cmap(cmap), norm=cNorm)

ax1.set_title(f'{Ref} error in Density from TTSE')
ax2.set_title(f'{Ref} error in Density from Bicubic')

for ax in [ax1, ax2]:

    ax.set_xlim(*T_lim)
    ax.set_ylim(*p_lim)

    ax.set_yscale('log')

    ax.grid(axis='both', which='both')
    ax.tick_params(axis='y', which='minor', left='off')

    ax.set_xlabel('Temperature [K]')
    ax.set_ylabel('Pressure [kPa]')

    ax.plot(Ts, ps, 'k', lw=4)

cbar_ax = fig.add_axes([0.85, 0.15, 0.06, 0.7])
CB = fig.colorbar(SC1, cax=cbar_ax)
CB.set_label(r'$(\rho/\rho_{EOS}-1)\times 100$ [%]')

plt.show()