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 https://github.com/CoolProp/CoolProp/issues/1301 https://github.com/CoolProp/CoolProp/issues/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

```python 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() ```