Multi-stage vapor injection impact on heat pump performance

[michox]

Hello dear community,

I am trying to calculate the impacts of a multi stage vapor injection on the performance of heat pumps, but unfortunately I have a hard time finding a configuration of boundary conditions that works. The goal is to limit the peak temperature (or peak pressure) and stay as close as possible to the 2-phase border during the compression to limit temperature and therefore energy consumption of the compression.

For N=0 the code works fine, but once I get to N=1 I keep getting a singularity in the jacobian matrix, and I can't figure out why, or which boundary conditions to choose to solve it.

Here is my code: https://github.com/Ecalia/heatpump-calculations/blob/main/VaporInjectionHeatPumpStudy.py

and the results can be seen / executed here https://github.com/Ecalia/heatpump-calculations/blob/main/heatpump-configurations.ipynb

Any help would be greatly appreciated!

Thanks in advance

[fwitte]

Hi @michox,

can you quickly add a flowchart to see, what the generic structure of the model will look like?

Also, I remember there was something with an injection compressor a long time ago (#230), maybe you can find some hints there already.

Best regards

[michox]

Hi @fwitte, thanks for the quick response, here is a flowchart for N=2:

I've seen the issue, but I was hoping to do it without specifying mass flow, which that solution seems to rely on?

[fwitte]

Ah okay, so you could (I'd recommend that actually) set the mass flow values and all those intermediate pressures for an initial guess, solve the network and then modify all the boundary conditions to you desired inputs and then resolve the system. That is a generally recommended approach, since the solver is very senstive to starting values. Once you have a stable starting point, you should be able to change whatever parameters you like :).

I drafted something, which is half ready, I'd provide you the second part tomorrow. It needs some modifications, i.e. upper and lower temperature values as user specifications and respecting critical point properties... (I have 100 °C as upper temperature, which would not work with R290).

from tespy.networks import Network
from tespy.components import (Compressor, CycleCloser, HeatExchangerSimple, Splitter, Merge, Valve)
from tespy.connections import Connection, Ref, Bus
from tespy.tools import document_models
from CoolProp.CoolProp import PropsSI


class InjectionHeatPump:

    def __init__(self, working_fluid, num_stages) -> None:
        self.working_fluid = working_fluid
        self.num_stages = num_stages
        self.setup_network()

    def setup_network(self):
        self.nwk = Network(fluids=[self.working_fluid], p_unit="bar", T_unit="C")

        condenser = HeatExchangerSimple("condenser")
        evaporator = HeatExchangerSimple("evaporator")
        cycle_closer = CycleCloser("cycle closer")

        splitters = {i: Splitter(f"splitter {i}") for i in range(self.num_stages - 1)}
        merges = {i: Merge(f"merge {i}") for i in range(self.num_stages - 1)}
        valves = {i: Valve(f"valve {i}") for i in range(self.num_stages)}
        compressors = {i: Compressor(f"compressor {i}") for i in range(self.num_stages)}

        c0 = Connection(condenser, "out1", cycle_closer, "in1", label="0")
        c1 = Connection(cycle_closer, "out1", valves[0], "in1", label="1")
        self.nwk.add_conns(c0, c1)

        for i in range(1, self.num_stages):
            _c1 = Connection(valves[i - 1], "out1", splitters[i - 1], "in1", label=f"vsp{i}")
            _c2 = Connection(splitters[i - 1], "out1", valves[i], "in1", label=f"spv{i}")
            _c3 = Connection(splitters[i - 1], "out2", merges[i - 1], "in2", label=f"spm{i}")
            self.nwk.add_conns(_c1, _c2, _c3)

        c2 = Connection(valves[len(valves) - 1], "out1", evaporator, "in1")
        c3 = Connection(evaporator, "out1", compressors[len(compressors) - 1], "in1")
        self.nwk.add_conns(c2, c3)

        for i in list(range(0, self.num_stages - 1))[::-1]:
            _c1 = Connection(compressors[i + 1], "out1", merges[i], "in1", label=f"cpm{i}")
            _c2 = Connection(merges[i], "out1", compressors[i], "in1", label=f"mcp{i}")
            self.nwk.add_conns(_c1, _c2)

        c4 = Connection(compressors[0], "out1", condenser, "in1")

        self.nwk.add_conns(c4)

        power_bus = Bus("power input")

        power_bus.add_comps(*[{"comp": compressor, "base": "bus"} for compressor in compressors.values()])
        self.nwk.add_busses(power_bus)

        [c.set_attr(eta_s=0.8) for c in compressors.values()]
        evaporator.set_attr(pr=1)
        condenser.set_attr(pr=1)

        p_high = PropsSI("P", "T", 100 + 273.15, "Q", 0, self.working_fluid) / 1e5
        c0.set_attr(fluid={self.working_fluid: 1}, x=0, p=p_high, m=2)
        p_low = PropsSI("P", "T", 0 + 273.15, "Q", 1, self.working_fluid) / 1e5
        c3.set_attr(x=1, p=p_low)
        print((p_low / p_high) ** (1 / self.num_stages))
        print(p_high)
        print(p_low)
        print(p_high * ((p_low / p_high) ** (1 / self.num_stages)) * ((p_low / p_high) ** (1 / self.num_stages)))

        self.nwk.get_conn("spv1").set_attr(p=2)
        self.nwk.get_conn("spv2").set_attr(p=.7)

        self.nwk.get_conn("spm1").set_attr(m=.2)
        self.nwk.get_conn("spm2").set_attr(m=.2)
        # self.nwk.get_conn("cpm1").set_attr(T=)

        self.nwk.solve("design")


InjectionHeatPump("R123", 3)

[fwitte]

If you are trying to build a generic model, this is some approach I would take for it:

from tespy.networks import Network
from tespy.components import (Compressor, CycleCloser, HeatExchangerSimple, Splitter, Merge, Valve)
from tespy.connections import Connection, Bus
import numpy as np
from CoolProp.CoolProp import PropsSI


class InjectionHeatPump:

    def __init__(self, working_fluid, num_stages) -> None:
        self.working_fluid = working_fluid
        self.num_stages = num_stages
        self.setup_network()

    def setup_network(self):
        self.nwk = Network(fluids=[self.working_fluid], p_unit="bar", T_unit="C")

        condenser = HeatExchangerSimple("condenser")
        evaporator = HeatExchangerSimple("evaporator")
        cycle_closer = CycleCloser("cycle closer")

        splitters = {i: Splitter(f"splitter {i}") for i in range(self.num_stages - 1)}
        merges = {i: Merge(f"merge {i}") for i in range(self.num_stages - 1)}
        valves = {i: Valve(f"valve {i}") for i in range(self.num_stages)}
        compressors = {i: Compressor(f"compressor {i}") for i in range(self.num_stages)}

        c0 = Connection(condenser, "out1", cycle_closer, "in1", label="0")
        c1 = Connection(cycle_closer, "out1", valves[0], "in1", label="1")
        self.nwk.add_conns(c0, c1)

        for i in range(1, self.num_stages):
            _c1 = Connection(valves[i - 1], "out1", splitters[i - 1], "in1", label=f"vsp{i}")
            _c2 = Connection(splitters[i - 1], "out1", valves[i], "in1", label=f"spv{i}")
            _c3 = Connection(splitters[i - 1], "out2", merges[i - 1], "in2", label=f"spm{i}")
            self.nwk.add_conns(_c1, _c2, _c3)

        c2 = Connection(valves[len(valves) - 1], "out1", evaporator, "in1")
        c3 = Connection(evaporator, "out1", compressors[len(compressors) - 1], "in1")
        self.nwk.add_conns(c2, c3)

        for i in list(range(0, self.num_stages - 1))[::-1]:
            _c1 = Connection(compressors[i + 1], "out1", merges[i], "in1", label=f"cpm{i}")
            _c2 = Connection(merges[i], "out1", compressors[i], "in1", label=f"mcp{i}")
            self.nwk.add_conns(_c1, _c2)

        c4 = Connection(compressors[0], "out1", condenser, "in1")

        self.nwk.add_conns(c4)

        power_bus = Bus("power input")

        power_bus.add_comps(*[{"comp": compressor, "base": "bus"} for compressor in compressors.values()])
        self.nwk.add_busses(power_bus)

        [c.set_attr(eta_s=0.8) for c in compressors.values()]
        evaporator.set_attr(pr=1)
        condenser.set_attr(pr=1)

        T_crit = PropsSI("Tcrit", self.working_fluid)
        p_high = PropsSI("P", "T", T_crit - 20, "Q", 0, self.working_fluid) / 1e5
        c0.set_attr(fluid={self.working_fluid: 1}, x=0, p=p_high, m=2)
        p_low = PropsSI("P", "T", T_crit - 100, "Q", 1, self.working_fluid) / 1e5
        c3.set_attr(x=1, p=p_low)
        # produces [p_high, ..., p_low]
        intermediate_pressures = np.geomspace(p_low, p_high, self.num_stages + 1)[::-1]

        for i in range(1, self.num_stages):
            # start: i = 1 -> intermediate_pressures[1] = first pressure below p_high
            self.nwk.get_conn(f"spm{i}").set_attr(p=intermediate_pressures[i])

        self.nwk.get_conn("spm1").set_attr(m=.2)
        self.nwk.get_conn("spm2").set_attr(m=.2)

        self.nwk.solve("design")

        self.nwk.print_results()


        for i in range(1, self.num_stages):
            # remove mass flow specifications again
            self.nwk.get_conn(f"spm{i}").set_attr(m=None)

        for i in range(0, self.num_stages - 1):
            # use compressor inlet state instead
            self.nwk.get_conn(f"mcp{i}").set_attr(x=1)
        self.nwk.solve("design")

        self.nwk.print_results()

InjectionHeatPump("R290", 3)

I put the compressor inlet states to saturated gas after the initial simulation. Of course, you could also remove the intermediate pressure values now and replace them by temperature. Or you can change upper and lower temperature bounds, ...

If you have more questions and some spare time, we can have a short chat via jit.si in the next days. Let me know, if you are interested.

Best

[michox]

Hey @fwitte thank you so much for this detailed response!

I have tried to implement the approach in my code and believe I have adjusted all the parameters equivalently and the first solving step runs just fine now, thanks!

Unfortunately, when I unset the massflow and specify the compressor inlet state instead and run it again, I get a RecursionError. Also, when I try to solve this for N=3 or above, I get a weird error that the starting composition of the fluid is empty before the first solution even runs...

Logs can be found here again https://github.com/Ecalia/heatpump-calculations/blob/main/heatpump-configurations.ipynb

I have taken your code and adjusted the labels so that it is easier to compare them with my code. The results of the first run are exactly the same.... only the cycle closer is at a different position, could that make a difference?

Also, is the bus used in any way here?

from tespy.networks import Network
from tespy.components import (Compressor, CycleCloser, HeatExchangerSimple, Splitter, Merge, Valve)
from tespy.connections import Connection, Bus
import numpy as np
from CoolProp.CoolProp import PropsSI


class InjectionHeatPump:

    def __init__(self, working_fluid, num_stages) -> None:
        self.working_fluid = working_fluid
        self.num_stages = num_stages
        self.setup_network()

    def setup_network(self):
        self.nwk = Network(fluids=[self.working_fluid], p_unit="bar", T_unit="C")

        condenser = HeatExchangerSimple("condenser")
        evaporator = HeatExchangerSimple("evaporator")
        cycle_closer = CycleCloser("cycle closer")

        splitters = {i: Splitter(f"splitter {i}") for i in range(self.num_stages - 1)}
        merges = {i: Merge(f"merge {i}") for i in range(self.num_stages - 1)}
        valves = {i: Valve(f"valve {i}") for i in range(self.num_stages)}
        compressors = {i: Compressor(f"compressor {i}") for i in range(self.num_stages)}

        c0 = Connection(condenser, "out1", cycle_closer, "in1", label="condenser-cycle_closer")
        c1 = Connection(cycle_closer, "out1", valves[0], "in1", label="cycle_closer-valve_1")
        self.nwk.add_conns(c0, c1)

        for i in range(1, self.num_stages):
            _c1 = Connection(valves[i - 1], "out1", splitters[i - 1], "in1", label=f"valve_{i}-splitter_{i}")
            _c2 = Connection(splitters[i - 1], "out1", valves[i], "in1", label=f"splitter_{i}-valve_{i+1}")
            _c3 = Connection(splitters[i - 1], "out2", merges[i - 1], "in2", label=f"splitter_{i}-merge_{i}")
            self.nwk.add_conns(_c1, _c2, _c3)

        c2 = Connection(valves[len(valves) - 1], "out1", evaporator, "in1", label=f"valve_{len(valves) }-evaporator")
        c3 = Connection(evaporator, "out1", compressors[len(compressors) - 1], "in1", label=f"evaporator-compressor_{len(compressors) }")
        self.nwk.add_conns(c2, c3)

        for i in list(range(0, self.num_stages - 1))[::-1]:
            _c1 = Connection(compressors[i + 1], "out1", merges[i], "in1", label=f"compressor_{i + 2}-merge_{i+1}")
            _c2 = Connection(merges[i], "out1", compressors[i], "in1", label=f"merge_{i+1}-compressor_{i+1}")
            self.nwk.add_conns(_c1, _c2)

        c4 = Connection(compressors[0], "out1", condenser, "in1", label=f"compressor_1-condenser")

        self.nwk.add_conns(c4)

        # power_bus = Bus("power input")

        # power_bus.add_comps(*[{"comp": compressor, "base": "bus"} for compressor in compressors.values()])
        # self.nwk.add_busses(power_bus)

        [c.set_attr(eta_s=0.8) for c in compressors.values()]
        evaporator.set_attr(pr=0.98)
        condenser.set_attr(pr=0.98)
        T_evap=-10
        T_cond=80
        p_high = PropsSI("P", "Q", 0, "T", 273.15 + T_cond, self.working_fluid) / 1e5
        p_low = PropsSI("P", "Q", 1, "T", 273.15 + T_evap, self.working_fluid) / 1e5


        c0.set_attr(fluid={self.working_fluid: 1}, x=0, p=p_high, m=2)
        c3.set_attr(x=1, p=p_low)
        # produces [p_high, ..., p_low]
        intermediate_pressures = np.geomspace(p_low, p_high, self.num_stages + 1)[::-1]

        for i in range(1, self.num_stages):
            # start: i = 1 -> intermediate_pressures[1] = first pressure below p_high
            print(f"setting pressure of splitter_{i}-merge_{i} to {intermediate_pressures[i]}")
            self.nwk.get_conn(f"splitter_{i}-merge_{i}").set_attr(p=intermediate_pressures[i], m=.2)


        self.nwk.solve("design")

        self.nwk.print_results()


        for i in range(1, self.num_stages):
            # remove mass flow specifications again
            self.nwk.get_conn(f"splitter_{i}-merge_{i}").set_attr(m=None)

        for i in range(0, self.num_stages - 1):
            # use compressor inlet state instead
            self.nwk.get_conn(f"merge_{i+1}-compressor_{i+1}").set_attr(x=1)
        self.nwk.solve("design")

        self.nwk.print_results()

InjectionHeatPump("R290", 2)

[fwitte]

The bus basically serves for easier calculation of the COP, it collects all compressors and their power input automatically. In the end you could have something like this:

hp = InjectionHeatPump("R290", 3)

print(hp.nwk.get_comp("condenser").Q.val / hp.nwk.busses["power input"].P.val)

On top of that, it is possible to add electrical/mechanical losses for busses.

To your code: It would take me a little bit longer going through that, because it is quite a lot already... I believe, you have two closed cycles (heat pump and the consumer?), that would mean you require two individual cycle closers (one per closed cycle). That might be the reason for the recursion error as well, but I am not sure, why it would only happen in the second solving step...

[michox]

Hi, thanks.

No, I am using a consumer cycle in other parts of my code, but not in this particular setup. Interestingly enough, even without any modification to the setup (i.e. commenting out the part that unsets the massflow) solving a second time results in that RecursionError... could that have anything to do with the connections being stored in a separate class member of my utility wrapper? I am not very familiar with pythons underlying magic... however, this only happens for this setup, my other heat pump setups solve without an issue when run multiple times...

[michox]

Got it to run! Turns out the issue was with the location of the cycle closer. After moving it to the same position where yours is, it works N=1 now, thanks! Is there information somewhere about why the position of the cycle closer matters?

It also will sometimes work for N=2 every other time, N=3 rarely and after running it like 20 times I also got it to run for N=4, I always get this error during the first run
WARNING:root:The starting value for the fluid composition of the connection merge_2-compressor_3 is empty. This might lead to issues in the initialisation and solving process as fluid property functions can not be called. Make sure you specified a fluid composition in all parts of the network.

whereas the connections for which this error message show up are not always the same, but there is always at least 4

Any ideas? it does not happen with your code...

[fwitte]

Great, the position of the cycle closer matters, people have been come across before also: #403. Although I am not sure 100 %, my guess is:

If you have a closed cycle, and specify mass flow anywhere, the cycle's mass flow will be overdetermined automatically due to the constraints of all components (also see https://tespy.readthedocs.io/en/main/basics/intro.html#establish-connections). If you have a combination of a merge and a splitter somewhere in the cycle, where either the merge's secondary mass flow (external mass flow going into the cycle) and the splitter's secondary mass flow (mass flow leaving the cycle) are not specified, those values will automatically equal out, effectively doing what the the cycle closer component does at this moment. That means, the cycle closer position should not be between the outlet of a merge and the inlet of a splitter, but between the outlet of a splitter and the inlet of a merge.

I do not have the time to test this out thouroughly, but I would be happy if you or the topic opener of #403 (@maltefritz) want to support a little bit here :).

To your second question:

In your notebook I found

WARNING:root:Reached maximum iteration count (50), calculation stopped. Residual value is 2.78e+02

That means, the network simulation has not converged, therefore the results are not meeting all constraints. There are several reasons, why this could happen. A frequent one is, if your specifications lead to states of the fluid, which are not within CoolProp's bounds for that fluid. That can lead to bad starting values for the next run and it even might modify your previous inputs, because some parameters are reprocessed wrongly in the simulation postprocessing... You might want to check on that.

Have a nice weekend!

[fwitte]

P.S.: There are two ways for you to know if the calculation has converged programatically:

• if nwk.lin_dep is True, that means there was a linear dependency and the solver stopped
• if nwk.res[-1] (res is a list storing the history of the residual per iteration) is larger than 1e-3, it means the calculation did not converge. Maybe adding a property nwk.converged to the network class might make sense here. I'll put that on my todo list.

[michox]

I figured out the problem with the
´WARNING:root:The starting value for the fluid composition of the connection merge_2-compressor_3 is empty. This might lead to issues in the initialization and solving process, as fluid property functions can not be called. Make sure you specified a fluid composition in all parts of the network.´
I had Water defined as part of the fluid network, but always set it to zero in my simplified cycles. Without it, my code now converges perfectly. Thanks @fwitte for the fast and thorough support, I really appreciate it, couldn't have done it without you!

[fwitte]

Good to hear that, one more tip if you still need water (or any other fluid e.g. for a different section of the cycle): provide starting values for the fluid at these points by specification of fluid0={...}. This will resolve that issue.

Thank you for reporting back, have a good weekend!