Help with connected heat exchangers modeling

[mortezaep]

Hi everyone

I am having difficulty modeling two connected heat exchangers. The fluids involved are water and R245fa. In the first heat exchanger, R245fa undergoes a phase change, and the code for this part works well, as shown below:

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

nw = Network(p_unit='bar', T_unit='C', h_unit='kJ / kg', iterinfo=False)

so1 = Source('source1')
si1 = Sink('sink1')
so2 = Source('source2')
si2 = Sink('sink2')
heatexchanger = HeatExchanger('heatexchanger')
heatexchanger.component()

so1_heatexchanger = Connection(so1, 'out1', heatexchanger, 'in1')
heatexchanger_si1 = Connection(heatexchanger, 'out1', si1, 'in1')
so2_heatexchanger = Connection(so2, 'out1', heatexchanger, 'in2')
heatexchanger_si2 = Connection(heatexchanger, 'out2', si2, 'in1')

nw.add_conns(so1_heatexchanger, heatexchanger_si1, so2_heatexchanger, heatexchanger_si2)

heatexchanger.set_attr(pr1=1, pr2=1, ttd_l=7)
so1_heatexchanger.set_attr(fluid={'water': 1}, T=140, p=40, m=70)

so2_heatexchanger.set_attr(fluid={'R245fa': 1}, x=0, p=3.8)
heatexchanger_si2.set_attr(x=1, p0=3.8)

# so2_heatexchanger.set_attr(fluid={'R245fa': 1}, x=0, m=135, p0=3.8)
# heatexchanger_si2.set_attr(x=1, p0=4.8)

nw.solve('design')
nw.print_results()

In the second code, I attempted to couple the first heat exchanger with another one, where R245fa remains in the liquid phase, and no phase change occurs. Unfortunately, the second code does not converge. The code for this scenario is provided below:

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

nw = Network(p_unit='bar', T_unit='C', h_unit='kJ / kg', iterinfo=False)

so1 = Source('source1')
si1 = Sink('sink1')
so2 = Source('source2')
si2 = Sink('sink2')
heatexchanger = HeatExchanger('heatexchanger')
heatexchanger2 = HeatExchanger('heatexchanger2')

heatexchanger.component()
heatexchanger2.component()

so1_heatexchanger = Connection(so1, 'out1', heatexchanger, 'in1')
heatexchanger_heatexchanger2 = Connection(heatexchanger, 'out1', heatexchanger2, 'in1')
heatexchanger2_si1 = Connection(heatexchanger2, 'out1', si1, 'in1')
so2_heatexchanger2 = Connection(so2, 'out1', heatexchanger2, 'in2')
heatexchanger2_heatexchanger = Connection(heatexchanger2, 'out2', heatexchanger, 'in2')
heatexchanger_si2 = Connection(heatexchanger, 'out2', si2, 'in1')

nw.add_conns(so1_heatexchanger, heatexchanger2_si1, so2_heatexchanger2, heatexchanger_si2, heatexchanger_heatexchanger2, heatexchanger2_heatexchanger)

heatexchanger.set_attr(pr1=1, pr2=1, ttd_l=7)
heatexchanger2.set_attr(pr1=1, pr2=1) #

so1_heatexchanger.set_attr(fluid={'water': 1}, T=140, p=40, m=70)
heatexchanger2_si1.set_attr(T=50, p0=40, m0=70) #

heatexchanger2_heatexchanger.set_attr(fluid={'R245fa': 1}, x=0, p=3.8, m0=130)
heatexchanger_si2.set_attr(x=1, p0=3.8, m0=130)

# heatexchanger2_heatexchanger.set_attr(fluid={'R245fa': 1}, x=0, m=135, p0=3.8)
# heatexchanger_si2.set_attr(x=1, p0=4.8)

nw.solve('design')
nw.print_results()

One workaround I found is to define two separate networks, which resolves the issue. This approach is also explained below:

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

nw = Network(p_unit='bar', T_unit='C', h_unit='kJ / kg', iterinfo=False)

so1 = Source('source1')
si1 = Sink('sink1')
so2 = Source('source2')
si2 = Sink('sink2')
heatexchanger = HeatExchanger('heatexchanger')
heatexchanger.component()

so1_heatexchanger = Connection(so1, 'out1', heatexchanger, 'in1')
heatexchanger_si1 = Connection(heatexchanger, 'out1', si1, 'in1')
so2_heatexchanger = Connection(so2, 'out1', heatexchanger, 'in2')
heatexchanger_si2 = Connection(heatexchanger, 'out2', si2, 'in1')

nw.add_conns(so1_heatexchanger, heatexchanger_si1, so2_heatexchanger, heatexchanger_si2)

heatexchanger.set_attr(pr1=1, pr2=1, ttd_l=7)
so1_heatexchanger.set_attr(fluid={'water': 1}, T=140, p=40, m=70)

so2_heatexchanger.set_attr(fluid={'R245fa': 1}, x=0, p=3.8)
heatexchanger_si2.set_attr(x=1, p0=3.8)

# so2_heatexchanger.set_attr(fluid={'R245fa': 1}, x=0, m=135, p0=3.8)
# heatexchanger_si2.set_attr(x=1, p0=4.8)

nw.solve('design')
nw.print_results()

nw2 = Network(p_unit='bar', T_unit='C', h_unit='kJ / kg', iterinfo=False)

so3 = Source('source3')
si3 = Sink('sink3')
so4 = Source('source4')
si4 = Sink('sink4')
heatexchanger2 = HeatExchanger('heatexchanger2')
heatexchanger2.component()

so3_heatexchange2 = Connection(so3, 'out1', heatexchanger2, 'in1')
heatexchanger2_si3 = Connection(heatexchanger2, 'out1', si3, 'in1')
so4_heatexchanger2 = Connection(so4, 'out1', heatexchanger2, 'in2')
heatexchanger2_si4 = Connection(heatexchanger2, 'out2', si4, 'in1')

nw2.add_conns(so3_heatexchange2, heatexchanger2_si3, so4_heatexchanger2, heatexchanger2_si4)

heatexchanger2.set_attr(pr1=1, pr2=1)

so3_heatexchange2.set_attr(fluid={'water': 1}, T=heatexchanger_si1.T.val, p=heatexchanger_si1.p.val, m=heatexchanger_si1.m.val)
heatexchanger2_si3.set_attr(T=50)

heatexchanger2_si4.set_attr(fluid={'water': 1}, T=so2_heatexchanger.T.val, p=so2_heatexchanger.p.val, m=so2_heatexchanger.m.val)

nw2.solve('design')
nw2.print_results()

Nonetheless, I believe there should be a way to fix the problem in the second code without resorting to separate networks. I would greatly appreciate any advice or suggestions to help me correct my mistake in designing the second code.

Kind regards
Morteza

[fwitte]

Hi Morteza,

thank you for reaching out. This is one of the issues you will often encounter while using tespy:

• You want to calculate the hot side first heat exchanger outlet temperature based on the evaporation temperature on the cold side.
• Since the cold side temperature is already fixed through p and x, it is only the T (and h) that should be determined on the hot side.
• The starting value for enthalpy (likely) is some value, that makes the state end up in two-phase region. Therefore, the derivative to enthalpy at the hot side connection between both heat exchangers will become zero, which it must not for the solver to find a solution.

There are two approaches to prevent this:

1. manually specify starting values in the correct locations
2. run two simulations in a row

The second approach is the one I usually take:

1. I build a model with as many (well educated guess) values for pressure and enthalpy (or p+T, p+x) directly specified and solve it once. Goal is, to just get a solution, that the solver was able to process. The solution itself does not need to be physically possible.
2. Unset all of those specifications, that were guesses and set the desired parameters instead.

If I want to go beyond scripting things, I will also create a class with a couple of methods, that I can control with an easy to use API (e.g. for workflow integration etc.). Unfortunately I have not yet found the time to create a tutorial on this topic, but it is on my list :D.

Solution with starting value

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

nw = Network(p_unit='bar', T_unit='C', h_unit='kJ / kg')

so1 = Source('source1')
si1 = Sink('sink1')
so2 = Source('source2')
si2 = Sink('sink2')
heatexchanger = HeatExchanger('heatexchanger')
heatexchanger2 = HeatExchanger('heatexchanger2')

heatexchanger.component()
heatexchanger2.component()

so1_heatexchanger = Connection(so1, 'out1', heatexchanger, 'in1')
heatexchanger_heatexchanger2 = Connection(heatexchanger, 'out1', heatexchanger2, 'in1')
heatexchanger2_si1 = Connection(heatexchanger2, 'out1', si1, 'in1')
so2_heatexchanger2 = Connection(so2, 'out1', heatexchanger2, 'in2')
heatexchanger2_heatexchanger = Connection(heatexchanger2, 'out2', heatexchanger, 'in2')
heatexchanger_si2 = Connection(heatexchanger, 'out2', si2, 'in1')

nw.add_conns(so1_heatexchanger, heatexchanger2_si1, so2_heatexchanger2, heatexchanger_si2, heatexchanger_heatexchanger2, heatexchanger2_heatexchanger)

heatexchanger.set_attr(pr1=1, pr2=1, ttd_l=7)
heatexchanger2.set_attr(pr1=1, pr2=1) #

so1_heatexchanger.set_attr(fluid={'water': 1}, T=140, p=40, m=70)
heatexchanger_heatexchanger2.set_attr(h0=250)
heatexchanger2_si1.set_attr(T=50)

heatexchanger2_heatexchanger.set_attr(fluid={'R245fa': 1}, x=0, p=3.8)
heatexchanger_si2.set_attr(x=1, p0=3.8)

nw.solve('design')

Two-step solution

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

nw = Network(p_unit='bar', T_unit='C', h_unit='kJ / kg')

so1 = Source('source1')
si1 = Sink('sink1')
so2 = Source('source2')
si2 = Sink('sink2')

heatexchanger = HeatExchanger('heatexchanger')
heatexchanger2 = HeatExchanger('heatexchanger2')

so1_heatexchanger = Connection(so1, 'out1', heatexchanger, 'in1')
heatexchanger_heatexchanger2 = Connection(heatexchanger, 'out1', heatexchanger2, 'in1')
heatexchanger2_si1 = Connection(heatexchanger2, 'out1', si1, 'in1')
so2_heatexchanger2 = Connection(so2, 'out1', heatexchanger2, 'in2')
heatexchanger2_heatexchanger = Connection(heatexchanger2, 'out2', heatexchanger, 'in2')
heatexchanger_si2 = Connection(heatexchanger, 'out2', si2, 'in1')

nw.add_conns(so1_heatexchanger, heatexchanger2_si1, so2_heatexchanger2, heatexchanger_si2, heatexchanger_heatexchanger2, heatexchanger2_heatexchanger)

so1_heatexchanger.set_attr(fluid={'water': 1}, T=140, p=40, m=70)
heatexchanger_heatexchanger2.set_attr(p=40, T=80)
heatexchanger2_si1.set_attr(p=40, T=50)

so2_heatexchanger2.set_attr(p=3.8)
heatexchanger2_heatexchanger.set_attr(fluid={'R245fa': 1}, x=0, p=3.8)
heatexchanger_si2.set_attr(x=1, p=3.8)

nw.solve('design')

heatexchanger.set_attr(ttd_l=7, pr1=1, pr2=1)
heatexchanger2.set_attr(pr1=1, pr2=1)

heatexchanger_heatexchanger2.set_attr(T=None, p=None)
heatexchanger2_si1.set_attr(p=None)
heatexchanger2_heatexchanger.set_attr(p=None)
heatexchanger_si2.set_attr(p=None)

nw.solve('design')

nw.print_results()

Possible solution for workflow integrations

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



class Model:

    def __init__(self):
        self._create_network()

    def _create_network(self):
        self.nw = Network(p_unit='bar', T_unit='C', h_unit='kJ / kg')

        so1 = Source('source1')
        si1 = Sink('sink1')
        so2 = Source('source2')
        si2 = Sink('sink2')

        heatexchanger = HeatExchanger('heatexchanger')
        heatexchanger2 = HeatExchanger('heatexchanger2')

        so1_heatexchanger = Connection(so1, 'out1', heatexchanger, 'in1', label='1')
        heatexchanger_heatexchanger2 = Connection(heatexchanger, 'out1', heatexchanger2, 'in1', label='2')
        heatexchanger2_si1 = Connection(heatexchanger2, 'out1', si1, 'in1', label='3')
        so2_heatexchanger2 = Connection(so2, 'out1', heatexchanger2, 'in2', label='11')
        heatexchanger2_heatexchanger = Connection(heatexchanger2, 'out2', heatexchanger, 'in2', label='12')
        heatexchanger_si2 = Connection(heatexchanger, 'out2', si2, 'in1', label='13')

        self.nw.add_conns(so1_heatexchanger, heatexchanger2_si1, so2_heatexchanger2, heatexchanger_si2, heatexchanger_heatexchanger2, heatexchanger2_heatexchanger)

        so1_heatexchanger.set_attr(fluid={'water': 1}, T=140, p=40, m=70)
        heatexchanger_heatexchanger2.set_attr(p=40, T=80)
        heatexchanger2_si1.set_attr(p=40, T=50)

        so2_heatexchanger2.set_attr(p=3.8)
        heatexchanger2_heatexchanger.set_attr(fluid={'R245fa': 1}, x=0, p=3.8)
        heatexchanger_si2.set_attr(x=1, p=3.8)

        self.nw.solve('design')

        heatexchanger.set_attr(ttd_l=7, pr1=1, pr2=1)
        heatexchanger2.set_attr(pr1=1, pr2=1)

        heatexchanger_heatexchanger2.set_attr(T=None, p=None)
        heatexchanger2_si1.set_attr(p=None)
        heatexchanger2_heatexchanger.set_attr(p=None)
        heatexchanger_si2.set_attr(p=None)

        self.nw.solve('design')
        self.stable_solution_path = 'stable_solution'
        self.nw.save(self.stable_solution_path)


    def get_param(self, obj, label, parameter):
        if obj == "Components":
            return self.nw.get_comp(label).get_attr(parameter).val
        elif obj == "Connections":
            return self.nw.get_conn(label).get_attr(parameter).val

    def set_params(self, **kwargs):
        if "Connections" in kwargs:
            for c, params in kwargs["Connections"].items():
                self.nw.get_conn(c).set_attr(**params)

        if "Components" in kwargs:
            for c, params in kwargs["Components"].items():
                self.nw.get_comp(c).set_attr(**params)

    def solve_model(self, **kwargs):
        self.set_params(**kwargs)

        try:
            self.nw.solve('design')
            if not self.nw.converged:
                self.nw.solve('design', init_only=True, init_path=self.stable_solution_path)
            else:
                # might need more checks here!
                if (
                        any(self.nw.results["HeatExchanger"]["Q"] > 0)
                    ):
                    self.solved = False
                else:
                    self.solved = True
        except ValueError as e:
            self.nw.lin_dep = True
            self.nw.reset_topology_reduction_specifications()
            self.nw.solve("design", init_only=True, init_path=self.stable_solution_path)

# %% Accessing and running the model

instance = Model()


specifications = {
    'Connections': {
        '1': {
            'T': 150
        }
    }
}

instance.solve_model(**specifications)
instance.nw.print_results()

specifications = {
    'Connections': {
        '1': {
            'T': None
        },
        '11': {
            'Td_bp': -20
        }
    }
}


instance.solve_model(**specifications)
instance.nw.print_results()