Allow to model mixtures of gases and liquids as well as CoolProp mixture back-end

[YoannUniKS]

Dear Francesco, dear all,

I'm using Tespy to calculate a heat exchanger and I'm running into a strange behaviour. I'm using version 0.6.3 of Tespy. Below a code example which should reproduce the problem.

The heat exchanger has water on the one side, and gas on the other side. I'm checking the result of the specific enthalpy calculation in the output "connections" csv file. When the gas is pure the specific enthalpy is consistent with the results of CoolProp. When the gas is a mixture, then the results differ.

As examples of the specific enthalpy calculations of the gas (secondary side) from the code below:

• pure methane: CoolProp and Tespy give the same results: 771,58 kJ/kg (5 °C) and 815,52 kJ/kg (20 °C)
• pure propane: CoolProp and Tespy give the same results: 216,99 kJ/kg (5 °C) and 254,23 kJ/kg (20 °C)
• mixture 50 % propane, 50 % methane: CoolProp gives 344,75 kJ/kg (5 °C) and 407,95 kJ/kg (20 °C) and Tespy gives 505,74 kJ/kg (5 °C) and 545,19 kJ/kg (20 °C)

I was wondering if I'm missing something or if the problem comes from Tespy.

I thank you in advance for the feedback.

from tespy.components import Sink, Source, HeatExchanger
from tespy.connections import Connection
from tespy.networks import Network

import CoolProp.CoolProp as CP

nw = Network(fluids=['water', 'methane', 'propane'], T_unit='C', p_unit='bar', m_unit='t / h', v_unit='m3 / h', h_unit='kJ / kg', iterinfo=False)

# Defining components
pri_in = Source('Primary side inlet')
pri_out = Sink('Primary side outlet')
sec_in = Source('Secondary side inlet')
sec_out = Sink('Secondary side outlet')
he = HeatExchanger('heat exchanger')

# Defining connections
pri_he = Connection(pri_in, 'out1', he, 'in1')
he_pri = Connection(he, 'out1', pri_out, 'in1')
sec_he = Connection(sec_in, 'out1', he, 'in2')
he_sec = Connection(he, 'out2', sec_out, 'in1')

nw.add_conns(pri_he, he_pri, sec_he, he_sec)

#Gas side
pu = 80 #bara
Tu = 5 #°C
Td = 20 #°C
dpg = 0 # bar

#Water side
pin = 2.5 #bara
Tin = 55 #°C
Tout = 45 #°C
dpw = 0 #bar

Vdot = 50

# set design attributes
he.set_attr(pr1=(pin-dpw)/pin, pr2=(pu-dpg)/pu, ttd_u=Tin-Td, design=['pr1', 'pr2', 'ttd_u'], offdesign=['zeta1', 'zeta2', 'kA_char'])

pri_he.set_attr(fluid={'water': 1, 'methane':0, 'propane':0}, p=pin, T=Tin, design=['v'])
sec_he.set_attr(fluid={'water': 0, 'methane':0.5, 'propane':0.5}, T=Tu, v=Vdot, p=pu)
he_pri.set_attr(T=Tout, design=['T'])

nw.solve(mode='design')
nw.save('design_test')

print('Design: \nVdot (gas): {a:0.1f} m³/h'.format(a=sec_he.v.val), '\nT (gas), in: {a:0.1f} °C, out: {b:0.1f} °C'.format(a=sec_he.T.val, b=he_sec.T.val), '\nVdot (water): {a:0.1f} m³/h'.format(a=pri_he.v.val), '\nT (water), in: {a:0.1f} °C, out: {b:0.1f} °C'.format(a=pri_he.T.val, b=he_pri.T.val), '\nkA: {a:0.1f} W/K'.format(a=he.kA.val))

fluid="methane[0.5]&propane[0.5]"
h_5 = CP.PropsSI('H','T', 5 + 273.15, 'P', 80 * 10**5, fluid)
h_20 = CP.PropsSI('H','T', 20 + 273.15, 'P', 80 * 10**5, fluid)

print('\nSpecific enthalpy gas inlet: {a:0.2f} kJ/kg \nSpecific enthalpy gas outlet: {b:0.2f} kJ/kg'.format(a=h_5/1000, b=h_20/1000))

[fwitte]

Hi,

enthalpy is a quantity with arbitrary reference point (there is no natural zero point). That means, if you want to compare enthalpies from different libraries, you should always use enthalpy differences subtracting a common reference point (temperature and pressure) from your enthalpy. Then these values should be roughly the same. TESPy uses simplfied ideal gas mixture rules for mixing different gases, CoolProp has more elborate routines, but they are quite fragile.

Best

[YoannUniKS]

Hi Francesco,

thank you for the quick reply. I agree with you, the difference is what matters. As the enthalpies are identical between CoolProp and Tespy for the pure components, I assumed that Tespy was also using CoolProp for mixtures, which is then not the case.

Back to the issue mentioned above. When comparing CoolProp and Tespy for a 50/50 propane/methane mixture defined above, the enthalpy difference is 63 kJ/kg with CoolProp and 39 kJ/kg with Tespy. That's far from a small difference (38 % less).

[fwitte]

Hmmm, that's unfortunate...

Two thoughts towards CoolProp

• "methane[0.5]&propane[0.5]" might refer to molar properties (not sure), but it does not explain the large difference.
• Not sure, whether to trust CoolProp mixture routines fully (e.g. https://github.com/CoolProp/CoolProp/issues/2333#issuecomment-1846169027). Usually, people point to using REFPROP in those instances. Do you happen to have a license and could check?

Appart from that: I think the issue is, that the saturation temperature of ethane at those pressures is very likely higher than the actual temperature. TESPy will account for only liquid or only gaseous phase, a condensation correction is currently only available for water (for performance reasons I do not want to check all mixture species by default). We can create a new mixture routine, which can account for condensation for all species of a mixture. Once I find the time for it, I will have a look and try to figure something out.

Slightly OT, but still relevant I think: We have validated mixture properties vs. Ebsilon Professional software, see https://github.com/KarimHShawky/Chemical-Exergy-in-TESPy/tree/main/validation/ebsilon. Those did match very good, also for higher pressure values at combustion, so I am wondering if they use the same property functions for mixtures as they are available through tespy. In that application condensation of water was relevant iirc.

Finally, I suggest you upgrade to the latest version of TESPy soon, because the implementation with the condensation check for generic mixtures will not happen in the older version, and the latest version also has a lot of other improvements, which you might benefit from :).

Have a nice weekend

[YoannUniKS]

Thanks for the feedback. Here some additional elements regarding your comments.

• it is indeed unclear in CoolProp if they use the molar or volume fraction. I did it the away round, considered 50/50 was the molar fraction (and the one required by CoolProp) and converted it to mass fraction for TESPy. It did not change much, I still get the same discrepancy.
• I do not have a REFPROP license, but I use to work with neqsim. Calculating the enthalpy difference with neqsim gives 62 kJ/kg, close to CoolProp (the extended code is below).
• I actually came to notice the problem by using TESPy to reproduce a real heat exchanger based on data sheet. While calculating the design point with TESPy I noticed that the heat transferred by the heat exchanger given in TESPy was lower than the one given in the data sheet (30 % less heat transferred calculated by TESPy compared to the data sheet). That is how I ended up comparing the enthalpy differences. So I think, this result combined with the two results above from neqsim and CoolProp tend to show that TESPy has an issue there.

Concerning your last remark regarding upgrading TESPy, I actually tried the latest version (0.7.1.post1), to check if maybe the problem was corrected there. But then I got an error on line 47, so I switched back to 0.6.3 (the error is linked to CoolProp). Here is the error message:

File "xxx.py", line 47, in nw.solve(mode='design') File "C:\Users\xxx\anaconda3\envs\tespy2\Lib\site-packages\tespy\networks\network.py", line 1963, in solve self.solve_loop(print_results=print_results) File "C:\Users\xxx\anaconda3\envs\tespy2\Lib\site-packages\tespy\networks\network.py", line 2017, in solve_loop self.solve_control() File "C:\Users\xxx\anaconda3\envs\tespy2\Lib\site-packages\tespy\networks\network.py", line 2293, in solve_control self.solve_components() File "C:\Users\xxx\anaconda3\envs\tespy2\Lib\site-packages\tespy\networks\network.py", line 2373, in solve_components cp.solve(self.increment_filter) File "C:\Users\xxx\anaconda3\envs\tespy2\Lib\site-packages\tespy\components\component.py", line 601, in solve self.residual[sum_eq:sum_eq + data.num_eq] = data.func( File "C:\Users\xxx\anaconda3\envs\tespy2\Lib\site-packages\tespy\components\heat_exchangers\base.py", line 585, in ttd_u_func T_o2 = o.calc_T(T0=o.T.val_SI) File "C:\Users\xxx\anaconda3\envs\tespy2\Lib\site-packages\tespy\connections\connection.py", line 742, in calc_T return T_mix_ph(self.p.val_SI, self.h.val_SI, self.fluid_data, self.mixing_rule, T0=T0) File "C:\Users\xxx\anaconda3\envs\tespy2\Lib\site-packages\tespy\tools\fluid_properties\functions.py", line 98, in T_mix_ph return inverse_temperature_mixture(kwargs) File "C:\Users\xxx\anaconda3\envs\tespy2\Lib\site-packages\tespy\tools\fluid_properties\helpers.py", line 77, in inverse_temperature_mixture return newton_with_kwargs( File "C:\Users\xxx\anaconda3\envs\tespy2\Lib\site-packages\tespy\tools\helpers.py", line 657, in newton_with_kwargs residual = target_value - function(function_kwargs) File "C:\Users\xxx\anaconda3\envs\tespy2\Lib\site-packages\tespy\tools\fluid_properties\mixtures.py", line 56, in h_mix_pT_ideal_cond return h_mix_pT_ideal(p, T, fluid_data, *kwargs) File "C:\Users\xxx\anaconda3\envs\tespy2\Lib\site-packages\tespy\tools\fluid_properties\mixtures.py", line 32, in h_mix_pT_ideal h += data["wrapper"].h_pT(pp, T) data["mass_fraction"] File "C:\Users\xxx\anaconda3\envs\tespy2\Lib\site-packages\tespy\tools\fluid_properties\wrappers.py", line 203, in h_pT self.AS.update(CP.PT_INPUTS, p, T) File "CoolProp\AbstractState.pyx", line 102, in CoolProp.CoolProp.AbstractState.update File "CoolProp\AbstractState.pyx", line 104, in CoolProp.CoolProp.AbstractState.update ValueError: For now, we don't support T [0 K] below Tmelt(p) [91.3316 K]

Code used for the calculations:

from tespy.components import Sink, Source, HeatExchanger
from tespy.connections import Connection
from tespy.networks import Network

from neqsim.thermo import fluid, fluid_df, TPflash, PSflash, PHflash

import CoolProp.CoolProp as CP

nw = Network(fluids=['water', 'methane', 'propane'], T_unit='C', p_unit='bar', m_unit='t / h', v_unit='m3 / h', h_unit='kJ / kg', iterinfo=False)

# Defining components
pri_in = Source('Primary side inlet')
pri_out = Sink('Primary side outlet')
sec_in = Source('Secondary side inlet')
sec_out = Sink('Secondary side outlet')
he = HeatExchanger('heat exchanger')

# Defining connections
pri_he = Connection(pri_in, 'out1', he, 'in1')
he_pri = Connection(he, 'out1', pri_out, 'in1')
sec_he = Connection(sec_in, 'out1', he, 'in2')
he_sec = Connection(he, 'out2', sec_out, 'in1')

nw.add_conns(pri_he, he_pri, sec_he, he_sec)

#Gas side
pu = 80 #bara
Tu = 5 #°C
Td = 20 #°C
dpg = 0 # bar

#Water side
pin = 2.5 #bara
Tin = 55 #°C
Tout = 45 #°C
dpw = 0 #bar

Vdot = 50

# set design attributes
he.set_attr(pr1=(pin-dpw)/pin, pr2=(pu-dpg)/pu, ttd_u=Tin-Td, design=['pr1', 'pr2', 'ttd_u'], offdesign=['zeta1', 'zeta2', 'kA_char'])

pri_he.set_attr(fluid={'water': 1, 'methane':0, 'propane':0}, p=pin, T=Tin, design=['v'])
sec_he.set_attr(fluid={'water': 0, 'methane':0.2668, 'propane':0.7332}, T=Tu, v=Vdot, p=pu)
he_pri.set_attr(T=Tout, design=['T'])

nw.solve(mode='design')
nw.save('design_test')

print('Design: \nVdot (gas): {a:0.1f} m³/h'.format(a=sec_he.v.val), '\nT (gas), in: {a:0.1f} °C, out: {b:0.1f} °C'.format(a=sec_he.T.val, b=he_sec.T.val), '\nVdot (water): {a:0.1f} m³/h'.format(a=pri_he.v.val), '\nT (water), in: {a:0.1f} °C, out: {b:0.1f} °C'.format(a=pri_he.T.val, b=he_pri.T.val), '\nkA: {a:0.1f} W/K'.format(a=he.kA.val))

fluid_CP="methane[0.5]&propane[0.5]"
h_5_CP = CP.PropsSI('H','T', 5 + 273.15, 'P', pu * 10**5, fluid_CP)
h_20_CP = CP.PropsSI('H','T', 20 + 273.15, 'P', pu * 10**5, fluid_CP)

print('CoolProp: \nEnthalpy difference: {a:0.2f} kJ/kg'.format(a=h_20_CP/1000 - h_5_CP/1000))

dic_fluid_NS = {'methane':0.5, 'propane':0.5}
fluid_NS = fluid("srk")
for key in dic_fluid_NS.keys():
    fluid_NS.addComponent(key, dic_fluid_NS[key], "mol/sec")

fluid_NS.setMixingRule("classic")
fluid_NS.setMultiPhaseCheck(True)

fluid_NS.setPressure(pu, "bara")

fluid_NS.setTemperature(5, "C")
TPflash(fluid_NS)
fluid_NS.initPhysicalProperties()
h_5_NS = fluid_NS.getEnthalpy("kJ/kg")

fluid_NS.setTemperature(20, "C")
TPflash(fluid_NS)
fluid_NS.initPhysicalProperties()
h_20_NS = fluid_NS.getEnthalpy("kJ/kg")

print('NeqSim:  \nEnthalpy difference: {a:0.2f} kJ/kg'.format(a=h_20_NS - h_5_NS))

Have a nice weekend too!

[fwitte]

Okay, thank you for that example. That seems to actually be a bug in the software, I found a fix but have to make sure, it does not interfere with other things. I'll open a pull request for the fix and I'll let you know, when it is available.

For the mixture properties: ethane is actually liquid at that pressure and that temperature, therefore the CoolProp and neqsim mxitures correctly take into consideration the partial evaporation when heating, tespy can only consider full or no evaporation in the mixture (the mixture properties are not built for liquid/gas mixtures for non pure fluids except water). I will also add a caution box to this section about it: https://tespy.readthedocs.io/en/main/modules/fluid_properties.html#tespy-fluid-properties-label.

With the latest version of tespy it might however be possible to directly integrate the CoolProp mixture backend into your models, I will create a separate pull request for that as well. Another idea could be to creat the custom functions for generic liquid/gas mixtures, but this might take a little bit more time. Finally, if you are interested, you could look into integrating neqsim into tespy. The new fluid property wrapper classes allow to integrate your own fluid property databases... You might also want to look into that, if you want to put a little bit of time and effort in :).

Thanks again for reporting!

Best

Francesco

[YoannUniKS]

Dear Francesco,

Thank you for the feedback. I appreciate the quick help!

I'm not quite sure why you mentioned ethane. In the mixture I defined there is only methane and propane.

But independently from that I agree, integrating neqsim into tespy would probably be the way. I looked into the page you linked, and it will require some time to get into the code of tespy. I was actually looking for a quick way to characterize my heat exchanger and that's how I came to tespy. As my heat exchanger is however not exactly the same as the one defined in tespy (mine is a mix of counter and cross flow, while in tespy it is a pure counter flow) I will invest the time I currently have in building a short model corresponding more precisely to my heat exchanger rather than in adding neqsim to tespy, for an approximate result in the end.

However in the medium term I will come back to your suggestion (integration of neqsim) and give it a try!

All the best!

[fwitte]

Hi @YoannUniKS, I misread that, thank you for clarifying!

If you are running more detailed calculations on heat exchangers, we could try to see, if those calculations could be integrated in tespy's heat exchanger models (or add new models repsectively). It is not that difficult to add new components actually :).

Have a nice weekend

[fwitte]

Issue with staring values for temperature (T0=0 Kelvin) resolved by https://github.com/oemof/tespy/pull/477