fwitte
Potential bug in Fluprodia's "calc_individual_isoline" method
Hello Francesco (@fwitte),
I want to report a possibly odd behavior of the Fluprodia package observed when creating the fluid property diagrams for a transcritical CO2 heat pump cycle modeled with TESPy. The heat pump cycle is the following:
Firstly, regarding the correctness and accuracy of my model, I believe there are no major issues, as the corresponding on-design simulation provides the expected results. Also, roughly the same results occur in the Streamlit dashboard of the "heatpumps" repository when similar heat pump parameter values are used (to the ones I use in the TESPy model) for this exact heat pump topology.
The issue is that when rerunning multiple times the same Python script (with the heat pump TESPy model), for several of these reruns an irregularity can be observed in the log(p)-h diagram of the refrigerant cycle (a redundant extended horizontal line of constant pressure appears in the diagram, as shown in the corresponding figure below). After diving deeper into this issue, I noticed that an extremely/abnormally high specific enthalpy value (at about 3458 kJ/kg, when the maximum h value for the refrigerant in the cycle is at about 517 kJ/kg) occurs at the compressor inlet (corresponds to the first value of the "h" np.array in the compressor datapoints dictionary produced by the Fluprodia's "calc_individual_isoline" method), at which location the fluid is in the saturated gas state. I have also noticed that this issue can only occur when any of the heat pump's state points (component inlet or outlet state) lie right above the saturation curve (saturated liquid or saturated gas curves). So, this could indicate an unstable behavior of the CoolProp TESPy backend in this region when using CO2 as the cycle's refrigerant.
However, the (on-design) TESPy results remain constant after each rerun. Additionally, the problem seems to appear only after the usage of Fluprodia's "calc_individual_isoline" method, which creates the "datapoints" dictionary (for the diagrams' isoline calculation). Considering these circumstances, I believe the issue may stem from a bug in the (corresponding) Fluprodia source code. If not, the problem should originate from the CoolProp library (so far, similar issues raised in the CoolProp forum have caught my attention, but none have been the same as this one, particularly concerning the CO2 fluid and the HEOS backend simultaneously).
Lastly, because this problem does not appear when rerunning the heat pump simulation through the environment of the "heatpumps" dashboard - which also uses TESPy and Fluprodia for the result generation - it is highly likely that the problem originates from my setup for the Fluprodia diagrams' generation. However, the code snippet in question seems to be correct.
Anyway, I have managed to circumvent the irregular property value issue in my code by using a monkey-patched version (which appears as a comment in the code below) of Fluprodia’s "calc_individual_isoline" method for automatically filtering out erroneous refrigerant property values (by defining appropriate bounds for all CO2 properties used in my diagrams). But of course, that technique does not confront the root cause of this problem.
I would be happy to provide any further explanation or disambiguation of the problem if needed.
Figures:
-
Normal
-
Faulty
Code:
from tespy.networks import Network
from tespy.components import (
HeatExchanger, CycleCloser, Compressor, Valve, Source, Sink
)
from tespy.connections import Connection, Bus
from CoolProp.CoolProp import PropsSI as PSI
working_fluid = "CO2" # cycle's refrigerant
nw = Network(T_unit="C", p_unit="bar", h_unit="kJ / kg")
# components
ev = HeatExchanger('evaporator')
gc = HeatExchanger('gas cooler')
cp = Compressor('centrifugal compressor')
va = Valve('expansion valve')
cc = CycleCloser('cycle closer')
sw_so = Source('source water source')
sw_si = Sink('source water sink')
dw_so = Source('district water source')
dw_si = Sink('district water sink')
# components' parameters
pr_dp_ev_hot = 0.98 # evaporator's hot stream design point pressure ratio (assuming its value)
pr_dp_ev_cold = 0.98 # evaporator's cold stream design point pressure ratio (assuming its value)
pr_dp_gc_hot = 0.98 # gas cooler hot stream design point pressure ratio (assuming its value)
pr_dp_gc_cold = 0.98 # gas cooler cold stream design point pressure ratio (assuming its value)
# connections
# source water loop
c11 = Connection(sw_so, "out1", ev, "in1", label="11") # in1 for heat exchanger's hot side inlet
c12 = Connection(ev, "out1", sw_si, "in1", label="12") # out1 for heat exchanger's hot side outlet
nw.add_conns(c11, c12)
# refrigerant (closed) loop
c0 = Connection(cc, "out1", gc, "in1", label="0")
c1 = Connection(gc, "out1", va, "in1", label="1")
c2 = Connection(va, "out1", ev, "in2", label="2")
c3 = Connection(ev, "out2", cp, "in1", label="3")
c4 = Connection(cp, "out1", cc, "in1", label="4")
nw.add_conns(c0, c1, c2, c3, c4)
# district water loop
c21 = Connection(dw_so, "out1", gc, "in2", label="21") # in2 for heat exchanger's cold side inlet
c22 = Connection(gc, "out2", dw_si, "in1", label="22") # out2 for heat exchanger's cold side outlet
nw.add_conns(c21, c22)
# connections' parameters
T_sw_so = 10 # source water supply temperature
T_dw_so = 40 # district heating network's return temperature
T_dw_si = 110 # district heating network's supply temperature
# power bus
motor = Bus("electric motor")
motor.add_comps(
{"comp": cp, "char": 1, "base": "bus"},
)
nw.add_busses(motor)
# parametrization
# making a good guess of pressure levels and enthalpy values
p_wf_GC_inlet = 184.69 # working fluid pressure in gas cooler inlet [bar], design (point) value,
# for CO2 it has to be P > Pcrit = 73.773 bar and T > Tcrit = 30.978 °C (for supercritical state)
T_wf_GC_outlet = 45 # working fluid temperature in gas cooler outlet [oC], assuming its value,
# for CO2 it has to be P > Pcrit = 73.773 bar and T > Tcrit = 30.978 °C (for supercritical state)
h_wf_GC_outlet = PSI('H', 'P', (pr_dp_gc_hot * p_wf_GC_inlet)*1e5, 'T', T_wf_GC_outlet + 273.15, working_fluid) # [J/kg]
#----------------------------------------------------------------------------------------------------------------------#
T_wf_ev = -4 # working fluid temperature in evaporator [oC], evaporation temperature, assuming its value,
# for CO2 it has to be < Tcrit = 30.978 °C (for subcritical state)
h_wf_ev_outlet = PSI('H', 'Q', 1, 'T', T_wf_ev + 273.15, working_fluid) # [J/kg], saturated gas (gas CO2)
p_wf_ev_outlet = PSI('P', 'Q', 1, 'T', T_wf_ev + 273.15, working_fluid) # [Pa], saturated gas (gas CO2)
#----------------------------------------------------------------------------------------------------------------------#
h_wf_ev_inlet = PSI('H', 'P', p_wf_ev_outlet / pr_dp_ev_cold, 'T', T_wf_ev + 273.15, working_fluid) # [J/kg], in case of two-phase region
# occurs in turbine outlet (temperature and pressure combination isn't enough for this thermodynamic state)
########################################################################################################################
gc.set_attr(Q=-48.4e6, pr1=pr_dp_gc_hot, pr2=pr_dp_gc_cold) # Q=48.4 MWth (Q=48.4e6), design (point) value (assuming pressure ratios)
ev.set_attr(pr1=pr_dp_ev_hot, pr2=pr_dp_ev_cold)
c1.set_attr(h=h_wf_GC_outlet / 1e3, fluid={working_fluid: 1}) # [kJ/kg]
#c2.set_attr(h=h_wf_ev_inlet / 1e3) # [kJ/kg], using enthalpy and pressure combination in case of two-phase region occurs
c3.set_attr(p=p_wf_ev_outlet / 1e5, h=h_wf_ev_outlet/ 1e3) # [bar], [kJ/kg]
c4.set_attr(T=141.65, p=p_wf_GC_inlet) # design (point) values
c11.set_attr(T=T_sw_so, p=1.013, fluid={'HEOS::Water': 1}) # design (point) values (assumption for pressure value)
c12.set_attr(T=7) # design (point) value
c21.set_attr(T=T_dw_so, p=10.204, fluid={'HEOS::Water': 1}) # design (point) values (assumption for pressure value)
c22.set_attr(T=T_dw_si) # design (point) value
print("====== 1st design run ======")
nw.solve("design")
nw.print_results()
print(f"Power input (electric) = {round(motor.P.val / 1e6, 3)} MWe")
print(f"Thermal output (heating capacity) = {round(abs(gc.Q.val)/1e6, 3)} MWth")
print(f"Thermal input (cooling capacity) = {round(abs(ev.Q.val)/1e6, 3)} MWth")
COP_hot = abs(gc.Q.val) / motor.P.val
COP_cold = abs(ev.Q.val) / motor.P.val
COP_total = COP_hot + COP_cold
print(f"Heating COP = {round(COP_hot, 3)}")
print(f"Cooling COP = {round(COP_cold, 3)}")
print(f"Total COP (excluding pumps' power, etc.) = {round(COP_total, 3)}", "\n")
# check
print(f"check result = {nw.lin_dep} | Target value: False", "\n") # if True the simulation run into a parametrisation's linear dependency
# unsetting stable starting values' specifications and setting target values instead - design point
c1.set_attr(h=None)
gc.set_attr(ttd_l=2) # assuming ... K [oC] pinch point (approach) temperature difference (ttd_l=...),
# assumption for the working fluid temperature in gas cooler outlet (common range: 35°C - 55°C),
# the approach temperature difference can be as low as 2 oC in compact HEXs (like PCHEs)
c2.set_attr(h=None, T=5) # assumption for the evaporation temperature (common range for CO2: -20°C - 20°C)
# pressure or temperature specification for two-phase region (for saturation curve / boiling line)
c3.set_attr(h=None, p=None, x=1) # assuming ... K [oC] of total/compressor superheat (superheated gas) at compressor suction
# (TESPy uses it as evaporator superheat),
# it depends on pressure corresponding to this connection
print("====== 2nd design run ======")
nw.solve("design")
nw.save("TCC HP_CO2_simple_design point")
nw.print_results()
print(f"compressor's eta_s = {cp.eta_s.val},", f"compressor's pr = {cp.pr.val},",
f"closed cycle's volumetric flow rate = {c1.v.val} (m^3)/s")
print(f"Power input (electric) = {round(motor.P.val / 1e6, 3)} MWe")
print(f"Thermal output (heating capacity) = {round(abs(gc.Q.val)/1e6, 3)} MWth")
print(f"Thermal input (cooling capacity) = {round(abs(ev.Q.val)/1e6, 3)} MWth")
COP_hot = abs(gc.Q.val) / motor.P.val
COP_cold = abs(ev.Q.val) / motor.P.val
COP_total = COP_hot + COP_cold
print(f"Heating COP = {round(COP_hot, 3)}")
print(f"Cooling COP = {round(COP_cold, 3)}")
print(f"Total COP (excluding pumps' power, etc.) = {round(COP_total, 3)}", "\n")
# check
print(f"check result = {nw.lin_dep} | Target value: False", "\n") # if True the simulation run into a parametrisation's linear dependency
# Plotting the heat pump's (refrigerant) on-design property diagrams using Fluprodia library
# Importing necessary library
import matplotlib.pyplot as plt
import numpy as np
from fluprodia import FluidPropertyDiagram
# Initial Setup
diagram = FluidPropertyDiagram("CO2")
diagram.set_unit_system(T='°C', p='bar', h='kJ/kg')
# # ******* monkey-patching (of Fluprodia’s "calc_individual_isoline" method): start *******
# # Saving original method
# original_calc_individual_isoline = FluidPropertyDiagram.calc_individual_isoline
#
# # Defining patch
# def safe_calc_individual_isoline(
# self,
# isoline_property=None,
# isoline_value=None,
# isoline_value_end=None,
# starting_point_property=None,
# ending_point_property=None,
# starting_point_value=None,
# ending_point_value=None
# ):
# datapoints = original_calc_individual_isoline(
# self,
# isoline_property=isoline_property,
# isoline_value=isoline_value,
# isoline_value_end=isoline_value_end,
# starting_point_property=starting_point_property,
# ending_point_property=ending_point_property,
# starting_point_value=starting_point_value,
# ending_point_value=ending_point_value
# )
#
# # Defining realistic physical bounds (for CO2 properties) for value filtering
# bounds = {
# 'h': (100, 700), # enthalpy in [kJ/kg] (range corresponding to the heat pump cycle specs)
# 'p': (10, 300), # pressure in [bar] (range corresponding to the heat pump cycle specs)
# 'T': (-50, 300), # temperature in [°C] (range corresponding to the heat pump cycle specs)
# 's': (500, 2500), # entropy in [J/kg·K] (range corresponding to the heat pump cycle specs)
# #'v': (1e-4, 5), # specific volume in [m³/kg] (range corresponding to the heat pump cycle specs)
# # -> not used in the diagrams / not sure about the values range
# }
#
# # Checking and sanitizing each property independently (filtering only out-of-bounds values, keeping NaNs untouched)
# for key, values in datapoints.items():
# if key not in bounds:
# continue
#
# low, high = bounds[key]
# values = np.array(values)
#
# # Identifying "nan" values (NaNs)
# nan_idx = np.isnan(values)
# if np.any(nan_idx):
# print(f"[Note] {key}: Found {nan_idx.sum()} NaN value(s): {values[nan_idx].tolist()}")
#
# # Identifying valid and invalid indices (keeping only the values being within bounds and NaN values)
# valid_idx = np.where(((values >= low) & (values <= high)) | nan_idx)[0]
# invalid_idx = np.where((values < low) | (values > high))[0]
#
# if len(valid_idx) < len(values):
# removed_vals = values[invalid_idx]
# print(f"[Filtered] {key}: Removed {len(removed_vals)} out-of-bounds value(s): {removed_vals.tolist()}")
#
# # Keeping all property arrays of the same heat pump component in sync length-wise
# for k in datapoints:
# if len(datapoints[k]) == len(values):
# datapoints[k] = np.array(datapoints[k])[valid_idx].tolist()
#
# return datapoints
#
# # Applying the patch
# FluidPropertyDiagram.calc_individual_isoline = safe_calc_individual_isoline
# # ******* monkey-patching: end *******
# Storing the heat pump's model result in the dictionary (generating plotting data)
result_dict = {}
# for every heat pump component in the closed (refrigerant) cycle
result_dict.update({ev.label: ev.get_plotting_data()[2]}) # because I need the data of evaporator's cold side stream
result_dict.update({cp.label: cp.get_plotting_data()[1]})
result_dict.update({gc.label: gc.get_plotting_data()[1]})
result_dict.update({va.label: va.get_plotting_data()[1]})
# Iterate over the results obtained from TESPy simulation
for key, data in result_dict.items():
# Calculate individual isolines for the desired diagrams
result_dict[key]['datapoints'] = diagram.calc_individual_isoline(**data)
print(f"{key} datapoints: {data['datapoints']}") # for debugging
# ***** for debugging - start *****
d = data['datapoints']
lens = {k: len(v) for k, v in d.items() if isinstance(v, (list, np.ndarray))}
print(f"{key} path lengths: {lens}")
# ***** for debugging - end *****
#============ Creating a figure and axis for plotting T-s diagram ============#
fig, ax = plt.subplots(1, figsize=(20, 10))
isolines = {
'Q': np.linspace(0, 1, 9),
'p': np.array([10, 20, 35, 50, 75, 100, 150, 200, 300, 500, 700]),
'v': np.array([]),
'h': np.arange(200, 650, 50)
}
# Set isolines for T-s diagram
diagram.set_isolines(**isolines)
diagram.calc_isolines()
# Draw isolines on the T-s diagram
diagram.draw_isolines(fig, ax, 'Ts', x_min=750, x_max=2500, y_min=-50, y_max=250)
# Adjust the font size of the isoline labels
for text in ax.texts:
text.set_fontsize(10)
# Plot T-s curves for each component
for key in result_dict.keys():
datapoints = result_dict[key]['datapoints']
_ = ax.plot(datapoints['s'], datapoints['T'], color='#ff0000', linewidth=2)
_ = ax.scatter(datapoints['s'][0], datapoints['T'][0], color='#ff0000')
# Set labels and title for the T-s diagram
ax.set_xlabel('Specific entropy, s in J/(kg*K)', fontsize=16)
ax.set_ylabel('Temperature, T in °C', fontsize=16)
ax.set_title('T-s Diagram of Heat Pump Cycle', fontsize=20)
# Set font size for the x-axis and y-axis ticks
ax.tick_params(axis='x', labelsize=12)
ax.tick_params(axis='y', labelsize=12)
plt.tight_layout()
# Save the T-s diagram plot as an SVG file
fig.savefig('TCC HP_CO2_simple_design_ts_diagram.svg')
#============ Creating a figure and axis for plotting log(p)-h diagram ============#
fig, ax = plt.subplots(1, figsize=(20, 10))
isolines = {
'Q': np.linspace(0, 1, 9),
'T': np.arange(-50, 250, 50),
'v': np.array([]),
's': np.arange(100, 3000, 200)
}
# Set isolines for log(p)-h diagram
diagram.set_isolines(**isolines)
diagram.calc_isolines()
# Draw isolines on the log(p)-h diagram
diagram.draw_isolines(fig, ax, 'logph', x_min=0, x_max=750, y_min=1e1, y_max=1e3)
# Adjust the font size of the isoline labels
for text in ax.texts:
text.set_fontsize(10)
# Plot log(p)-h curves for each component
for key in result_dict.keys():
datapoints = result_dict[key]['datapoints']
_ = ax.plot(datapoints['h'], datapoints['p'], color='#ff0000', linewidth=2)
_ = ax.scatter(datapoints['h'][0], datapoints['p'][0], color='#ff0000')
#print(f"{key} h: {datapoints['h']}") # for debugging
#print(f"{key} p: {datapoints['p']}") # for debugging
# Set labels and title for the log(p)-h diagram
ax.set_xlabel('Specific enthalpy, h in kJ/kg', fontsize=16)
ax.set_ylabel('Pressure, p in bar', fontsize=16)
ax.set_title('log(p)-h Diagram of Heat Pump Cycle', fontsize=20)
# Set font size for the x-axis and y-axis ticks
ax.tick_params(axis='x', labelsize=12)
ax.tick_params(axis='y', labelsize=12)
plt.tight_layout()
# Save the log(p)-h diagram plot as an SVG file
fig.savefig('TCC HP_CO2_simple_design_logph_diagram.svg')
User specifications OS: Windows 10 Home Python: 3.11 Python IDE: PyCharm Community Edition 2024.3.5 TESPy: 0.7.9 CoolProp: 6.7.0 Fluprodia: 3.5.1
Hi @MikeSennis,
this is indeed a bug. I had noticed it from time to time already, it seems to appear, when the first point is at saturation or just above that. A fix could be to clean the data for these kind of large jumps. What you can do, is just remove the outlier manually in the .plot() function call.
Best
