Thermal Storage as Component

[jubra97]

Hello,
thank you for your helpful package. I’m currently working on a model of a network to distribute heat between an ORC power plant and a district heating network. In this network a thermocline thermal storage is included as well. But I have a hard time to model this storage with your framework.

My first goal is to model the storage in a way that the outgoing mass flows have a fixed temperature if the storage in not full or empty. The ingoing temperatures can be ignored. There must be no flow in both incoming mass flow streams at the same time. One of the incoming flows must be zero. The pressure at the bottom of the storage tank must be 1.8 bar higher than the pressure at the top. To compute the current hot water level in the tank I use multiple simulations and integrate the inflows and outflows. This configuration is shown in the figure below. The green color at the connections indicate that these values are the model constrains. The black values should be calculated.

And here is the corresponding code:

from tespy.components.component import Component
from tespy.tools.data_containers import ComponentProperties as dc_cp
from tespy.networks import Network
from tespy.components import Source, Sink
from tespy.connections import Connection
from CoolProp.CoolProp import PropsSI
import numpy as np

class TES(Component):
    def __init__(self, label, step_time=5, **kwargs):  # time in minutes
        self.step_time = step_time
        self.hot_site_temp = 160
        self.cold_site_temp = 60
        super().__init__(label, **kwargs)

    @staticmethod
    def component():
        return 'TES'

    @staticmethod
    def inlets():
        return ["in1", "in2"]

    @staticmethod
    def outlets():
        return ["out1", "out2"]

    def get_variables(self):
        return {
            'hot_water_volume': dc_cp(min_val=0),
            'tank_volume': dc_cp(min_val=0),
            'hot_water_level': dc_cp(min_val=0, max_val=1)

        }

    def get_mandatory_constraints(self):
        return {
            'mass_flow_constraints': {
                'func': self.mass_flow_func, 'deriv': self.mass_flow_deriv,
                'constant_deriv': True, 'latex': self.mass_flow_func_doc,
                'num_eq': self.num_i},
            'fluid_constraints': {
                'func': self.fluid_func, 'deriv': self.fluid_deriv,
                'constant_deriv': True, 'latex': self.fluid_func_doc,
                'num_eq': self.num_nw_fluids * self.num_i},
            'thermal_constraints': {
                'func': self.thermal_func,
                'deriv': self.thermal_deriv,
                'constant_deriv': False, 'latex': self.thermal_func_doc,
                'num_eq': 2},
        }

    def thermal_func(self):
        r"""
        Equation for thermal output calculation.
        """
        res = []

        res += [(self.outl[0].h.val_SI - PropsSI("H", "T", self.hot_site_temp + 273.15, "P", self.outl[0].p.val_SI, "Water"))]
        res += [(self.outl[1].h.val_SI - PropsSI("H", "T", self.cold_site_temp + 273.15, "P", self.outl[1].p.val_SI, "Water"))]
        return np.array(res)

    def thermal_func_doc(self, label):
        return ""

    def thermal_deriv(self, increment_filter, k):
        r"""
        Calculate partial derivatives of thermal output.

        Parameters
        ----------
        increment_filter : ndarray
            Matrix for filtering non-changing variables.

        k : int
            Position of derivatives in Jacobian matrix (k-th equation).
        """
        if not increment_filter[2, 2]:
            self.jacobian[k, 2, 2] = self.numeric_deriv(self.thermal_func, 'h', 2)[0]
            
        if not increment_filter[3, 2]:
            self.jacobian[k, 3, 2] = self.numeric_deriv(self.thermal_func, 'h', 3)[1]

 
    def calc_parameters(self):
        r"""Postprocessing parameter calculation."""
        if abs(self.inl[0].m.val_SI) > 1e-3 and abs(self.inl[1].m.val_SI) > 1e-3:
            raise ValueError("Flow is in only one inlet allowed.")
        self.hot_water_volume.val = self.hot_water_volume.val + (
                    (self.inl[1].m.val_SI * (self.step_time * 60)) / PropsSI("D", "T", self.hot_site_temp + 273.15, "P",
                                                                             self.inl[1].p.val_SI, "Water"))
        self.hot_water_level.val = self.hot_water_volume.val / self.tank_volume.val


if __name__ == "__main__":
    fluid_list = ["Water"]
    network = Network(fluids=fluid_list)
    network.set_attr(p_unit="bar", T_unit="C", h_unit="kJ / kg", m_unit="kg / s", m_range=[-0.1, 20])

    src_cold = Source("Src Cold")
    src_hot = Source("Src Hot")
    sink_cold = Sink("Sink Cold")
    sink_hot = Sink("Sink Hot")
    tes = TES("TES")

    c1 = Connection(src_cold, "out1", tes, "in1", label="c1")
    c2 = Connection(tes, "out1", sink_hot, "in1", label="c2")
    c3 = Connection(src_hot, "out1", tes, "in2", label="c3")
    c4 = Connection(tes, "out2", sink_cold, "in1", label="c4")

    network.add_conns(c1, c2, c3, c4)
    tes.set_attr(hot_water_volume=0.5, tank_volume=np.pi * 2 ** 2 * 10, hot_water_level=0)
    c1.set_attr(fluid={"Water": 1}, p=9.3, m=0, T=50)
    c3.set_attr(fluid={"Water": 1}, p=8, m=3, T=140)
    c2.set_attr(p=8)
    c4.set_attr(p=9.3)

    network.solve("design", max_iter=200, init_only=False)

I tried to use the HeatExchanger module as a guide. My two main problems are with the thermal_deriv function, where I’m not sure if I’ve placed the derivatives in the right place in the Jacobian matrix and with calculation the output pressure. If I try to run my code, I get a singularity in the Jacobian matrix. If I use the pseudo inverse instead of the normal inverse my code runs but sometimes produces some strange bugs. Could you help me with the required equations?
Thank you in advance!

[fwitte]

Dear @jubra97,

thank you for your idea and for reaching out! I have been on holidays, and did not yet have a lot of time to look into your suggestion, I will hopefully be able to do that on the next weekend. Quick general comments after scimming through your comment:

• Massflows equal to zero cannot be modelled well in TESPy, therefore I'd suggest to
• separate the charging and discharging functionalities for the TESPy side, i.e. creating two separate TESPy models, one for the storage charging mode and one for the discharging mode

I'll come back to you later

Best

[jubra97]

Hello fwitte,
Thank you for you advice! I think I found a way to solve the problem without using two separated models. To do this I created a subsystem with two HeatExchangerSimple and some extra logic. This logic is responsible for controlling the current level of the storage and to prevent that the storage is in charge and discharge mode at the same time. This basically means that there is only mass flow allowed in one direction. To control the output temperatures of the two HeatExchangerSimple I set the temperature at the output connection manually. After every simulation I call the extra logic in the subsystem to calculate the level of the tes and to check if anything is out of bounds. I attached an image of the model and the corresponding code:

from tespy.components import Merge, Splitter, Subsystem, HeatExchangerSimple
from tespy.connections import Connection
import numpy as np
from CoolProp.CoolProp import PropsSI

class TESSystem(Subsystem):
    """
    Class that represents the TES.
    """

    def __init__(self, label: str, step_time: float=5, tank_diameter: float = 3, tank_height: float = 10, tank_max_change_vel: float = 0.001, hot_water_volume: float = None):
        super().__init__(label)
        self.step_time = step_time  # time in minues
        self.tank_diameter = tank_diameter
        self.tank_height = tank_height
        self.tank_volume = np.pi * (self.tank_diameter / 2) ** 2 * self.tank_height
        self.tank_max_change_vel = tank_max_change_vel
        self.max_flow = (np.pi * (self.tank_diameter / 2) ** 2 * self.tank_max_change_vel) * 1000  # approx for density of water
        self.hot_water_volume = hot_water_volume if hot_water_volume else self.tank_volume / 2
        self.hot_water_level = self.hot_water_volume / self.tank_volume
        self.set_default_component_constraints()
        
    def create_comps(self):
        self.comps["hx_charge"] = HeatExchangerSimple("TES Charge")
        self.comps["hx_discharge"] = HeatExchangerSimple("TES Discharge")

    def check_constrains(self):
        if abs(self.comps["hx_charge"].inl[0].get_attr("m").val_SI) > 1e-3 and abs(self.comps["hx_discharge"].inl[0].get_attr("m").val_SI) > 1e-3:
            raise ValueError("Only massflow in one direction allowed.")
        if self.comps["hx_charge"].inl[0].get_attr("m").val_SI > self.max_flow or self.comps["hx_discharge"].inl[0].get_attr("m").val_SI > self.max_flow:
            raise ValueError("Charge velocity is too high.")
        if self.hot_water_level > 1:
            raise ValueError("TES is already full!")
        if self.hot_water_level < 0:
            raise ValueError("TES is already empty!")

    def update_storage(self):
        self.hot_water_volume = self.hot_water_volume + (
                    ((self.comps["hx_charge"].inl[0].m.val_SI - self.comps["hx_discharge"].inl[0].m.val_SI) * (self.step_time * 60)) / PropsSI("D", "T", 160 + 273.15, "P", self.comps["hx_charge"].inl[0].p.val_SI, "Water"))
        self.hot_water_level = self.hot_water_volume / self.tank_volume

    def postprocess(self):
        self.update_storage()
        self.check_constrains()

    def set_default_component_constraints(self):
        self.comps["hx_charge"].get_attr("pr").max_val = 2
        self.comps["hx_charge"].get_attr("zeta").min_val = -self.comps["hx_charge"].get_attr("zeta").max_val

        
if __name__ == "__main__":
    from tespy.networks import Network
    from tespy.components import Sink, Source, HeatExchangerSimple
    from tespy.connections import Ref
    fluid_list = ["Water"]
    network = Network(fluids=fluid_list)
    network.set_attr(p_unit="bar", T_unit="C", h_unit="kJ / kg", m_unit="kg / s", m_range=[-0.1, 20])

    he = HeatExchangerSimple("HE")
    tes = TESSystem("TES")
    merge_j = Merge("Merge J")
    splitter_j = Splitter("Splitter J")
    merge_g = Merge("Merge G")
    splitter_e = Splitter("Splitter E")
    sink = Sink("Sink")
    source = Source("source")

    c_27 = Connection(merge_g, "out1", splitter_e, "in1")
    c_16 = Connection(merge_j, "out1", sink, "in1")
    c_4b = Connection(source, "out1", merge_g, "in2")
    c_15_charge = Connection(splitter_j, "out1", tes.comps["hx_charge"], "in1")
    c_15_intern = Connection(splitter_j, "out2", merge_j, "in1")
    c_15_discharge = Connection(tes.comps["hx_discharge"], "out1", merge_j, "in2")
    c_14_charge = Connection(tes.comps["hx_charge"], "out1", merge_g, "in1")
    c_14_discharge = Connection(splitter_e, "out2", tes.comps["hx_discharge"], "in1")
    c_10 = Connection(splitter_e, "out1", he, "in1")
    c_6 = Connection(he, "out1", splitter_j, "in1")

    network.add_subsys(tes)
    network.add_conns(c_27, c_16, c_4b, c_15_charge, c_15_intern, c_15_discharge, c_14_charge, c_14_discharge, c_10, c_6)
    
    he.set_attr(Q=3.5e6)

    c_27.set_attr(fluid={"Water": 1})
    c_15_discharge.set_attr(T=160)
    c_14_charge.set_attr(T=60)
    c_14_discharge.set_attr(p0=9.3)
    c_14_charge.set_attr(m=0, p=Ref(c_15_discharge, 1, tes.tank_height/10))
    c_15_discharge.set_attr(m=1, p=8)
    c_4b.set_attr(T=58)
    c_10.set_attr(m=8)

    network.solve("design")
    network.print_results()
    tes.update_storage()

Come back to me if you have any questions to the code.

It would by nice to get this subsystem working as a component. Dou you have any idea how it would be possible to transfer this logic in a component?

Best

[fwitte]

@jubra97,

thank you for updating on this. I am struggling to find the time to go deep into this. Generally, I would think, that this is possible, but in some way it is not the core of what TESPy does. Nevertheless, coming up with some kind of template to reuse your structure/logic for similar components has some value for other people I think. Maybe we can meet some time in the future and discuss this? There is the TESPy user meeting, every 3rd Monday of the month at 17:00 CEST, if you want you can join for the August meeting.

Best

[jubra97]

Sorry for the late answer. Would that be today?
Regards

[fwitte]

That is today, yes!

[jgunstone]

I was searching for some representation of thermal storage in TESpy and came across this -
I'd like to model heat-pump operation for large-scale commercial building (or central plant for residential blocks).
A key component of such systems is thermal storage to prevent cycling of the heat pump.
from my POV it would be great if there was a built in component to handle this - or else some documentation describing a suggested approach.

many thanks

[fwitte]

Hi @jgunstone,

thank you for reaching out. There was actually an approach recently implemented by someone in the discussions (#505, title does not show that immediately...).

I would like to put this on the list for more examples/tutorials in the docs (see #506). Maybe we can discuss the concept for a simple to replicate and understand integration of storage into a tespy simulation model in on of the next online user meeting, would you like to join?

Best

Francesco

[jgunstone]

hi @fwitte - I'll take a look at #505. Yes I'd be very interested to join the next meeting - thanks

John