Issue with TTSE interpolation when using setState_pT for supercritical CO2

[casella]

I just tried out the newly released 3.3.1 version using Dymola 2022 under Windows 10, 64 bits. The examples in the ExternalMedia library seem to work fine.

I then tested the BICUBIC and TTSE CO2 medium model. I confirm that the issue I had in #29 is now fixed, the first time I ran the model it took some time to generate the tables, and then the simulation ran correctly. Subsequent simulations went like a breeze.

I tested the CO2 medium model on an isobaric line at 90 bars, which is supercritical, but not too much, so I expect to see significant real-gas effects. The setState_pT function is used, see the test model TestInterpolatedCO2.mo.txt. I report the plots of density:

There is clearly a problem with the TTSE interpolation between 85 and 150 bars. Above 150 bars, all seems fine. @jowr, could you please have a look?

[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()