Open loop DH network with heating consumer components

[mike280512]

Hi all,
I'm in the process of modelling an existing DH network which is an open loop system using a geothermal resource. Heated water is extracted in a source location and then returned at another sink location. A network of supply pipes transports heated water from the source, through customer heat exchangers, before return pipes collect returned water to be injected at the sink location.

It is important to model both supply and return pipes, since they are not identical layouts and have different junction heights.

I would like to use pandapipes to model this system, being able to control the known mass flows through each customer heat exchanger. However, it seems that, unless the network is a closed loop, it is not possible to use a mass flow controlling component, e.g. a heat consumer, in an open loop system - the pipeflow always fails to converge.

As a minimum working example, the following will fail to converge:

net = pp.create_empty_network(fluid="water")
net.ambient_temperature_k = 293.15 

junction_1 = pp.create_junction(net, 1, 298, 0, "junction_1")
junction_2 = pp.create_junction(net, 1, 298, 0, "junction_2")
junction_3 = pp.create_junction(net, 1, 298, 0, "junction_3")

source = pp.create_ext_grid(net,junction_1,p_bar=1,t_k=298)
pipe = pp.create_pipe_from_parameters(net, junction_1, junction_2, diameter_m=0.1, length_km=0.1, k_mm=0.2, loss_coefficient=0)
heat_consumer = pp.create_heat_consumer(net, junction_2, junction_3, qext_w=-500, controlled_mdot_kg_per_s=0.5)
sink = pp.create_ext_grid(net,junction_3,t_k=298)

pp.pipeflow(net)

I've tried different combinations of component to represent the sink, including an external grid (as above), a sink with matching mass flow rate, or simply a junction, all with no success.

Splitting the supply and return into two separate networks works well for the supply side, but is a less elegant solution than I hoped would be possible:

net = pp.create_empty_network(fluid="water")
net.ambient_temperature_k = 293.15  # Set ambient temperature

junction_1 = pp.create_junction(net, 1, 298, 0, "junction_1")
junction_2 = pp.create_junction(net, 1, 298, 0, "junction_2")
junction_3 = pp.create_junction(net, 1, 298, 0, "junction_3")

source = pp.create_ext_grid(net,junction_1,p_bar=1,t_k=298)
pipe = pp.create_pipe_from_parameters(net, junction_1, junction_2, diameter_m=0.1, length_km=0.1, k_mm=0.2, loss_coefficient=0)
heat_exchanger = pp.create_heat_exchanger(net, junction_2, junction_3, qext_w=-500, diameter=0.02, loss_coefficient=1)
to_return = pp.create_sink(net,junction_3,mdot_kg_per_s=0.5)
pp.pipeflow(net)

return_net = pp.create_empty_network(fluid="water")
return_net.ambient_temperature_k = 293.15  # Set ambient temperature

junction_1 = pp.create_junction(return_net, 1, 298, 0, "junction_1")
junction_2 = pp.create_junction(return_net, 1, 298, 0, "junction_2")

return_pressure = net.res_heat_exchanger['p_to_bar'][0]
return_temp = net.res_heat_exchanger['t_to_k'][0]

from_supply = pp.create_source(return_net,junction_1,mdot_kg_per_s=0.5)
pipe_return = pp.create_pipe_from_parameters(return_net, junction_1, junction_2, diameter_m=0.1, length_km=0.1, k_mm=0.2, loss_coefficient=0)
sink = pp.create_ext_grid(return_net,junction_2,p_bar=1)

pp.pipeflow(return_net)

Can the developers offer any advice on this? Is this the expected behaviour, i.e. mass flow can only be controlled at the sinks for an open loop network?

I can 'artificially' connect the source and sink with a constant pressure pump to achieve convergence in a closed loop, but I am interested in determining what the pressure at the sink is. Supplying both constant pressure and pressure lift prevents this pressure from being calculated.

[mike280512]

I think I've resolved this now - the problem above is under-constrained when using the heat consumer component, the pressure downstream must be fixed in some way. Specifying the pressure for the 'sink' external grid returns a converged solution.

The end result is a bit like the separate supply and return network example, in that the pressure is fixed at the sink external grid and determines the pressure for the 'to junction' of the heat consumer. I had expected that the pressure profile would be calculated from source to sink, in the direction of mass flow, but clearly this is not the case.

[jnettels]

Hi there, just a random user here. I will probably not directly help you solve your issue, but I wanted to contribute my own thoughts here (since my problem seems very similar) instead of opening a new topic.
You said you are "interested in determining what the pressure at the sink is". While I am trying to simulate a closed loop system, I guess it boils down to the same question: If a pump defines the forward flow pressure and the pressure drops along the pipes and at each consumer substation, what is the return pressure at the pump?

My code example grew a little bit less minimal but allows to choose a variable number of consumers and combine "heat_exchanger" plus "flow_control" plus "valve" instead of "heat_consumer" components. The pump defines a mass flow, as well as all consumers except one to avoid a singular matrix. With pandapipes 0.13 'bidirectional' mode worked, but now in 0.14 only 'hydraulics' works, which might be another issue. Anyway, this is the net with 4 consumers:

The next plot shows the resulting pressure at each junction. My red arrows show my perception about how the pressure information "travels" in this system. The pressure in the lower (return) pipe segment is determined by the first connection to the forward segment. Then pressure increases along the return segment (against the flow direction) due to friction losses. But we see that the return pressure after a consumer can become higher then the forward pressure, which seems very problematic.

The problem can be summarized in this plot of pressure vs. distance, where the return pressure should not be higher than the forward/supply pressure.

If using closed-loop systems is encouraged, am I using the library incorrectly?

This is my code:

Spoiler warning

import pandas as pd
import pandapipes as pp
from matplotlib import pyplot as plt

# create empty net
net = pp.create_empty_network(fluid ="water")

pn_bar=8
t_flow_k = 273.15+70
cp = 4.18 # [kJ/(kg K)]
dT = 30  # K
# pressure_loss_coefficient = 10000
pressure_loss_coefficient = 0
n_consumers = 4

# Define consumers with required thermal power and mass flow
consumers = pd.DataFrame(dict(qext_w=[20000]*n_consumers))
consumers['mdot_kg_per_s'] = consumers['qext_w'] / (cp * dT) / 1000

# Create pump and its junctions
j0 = pp.create_junction(net, pn_bar=pn_bar, tfluid_k=293.15, name="junction 0", geodata=(0, 0))
j1 = pp.create_junction(net, pn_bar=pn_bar, tfluid_k=293.15, name="junction 1", geodata=(0, 30))

pp.create_circ_pump_const_mass_flow(
    net, return_junction=j0, flow_junction=j1, p_flow_bar=pn_bar,
    mdot_flow_kg_per_s=consumers['mdot_kg_per_s'].sum(),
    t_flow_k=t_flow_k,
    type="auto", name="const_mass_pump", index=None, in_service=True)

# Test external_grid plus sink, instead of pump
# ext_grid = pp.create_ext_grid(net, junction=j1, p_bar=pn_bar, t_k=t_flow_k, name="Grid Connection")
# sink = pp.create_sink(net, junction=j0, mdot_kg_per_s=consumers['mdot_kg_per_s'].sum(), name="Sink")

# Create a variable number of consumers and their junctions
js_f = [j1]
js_r = [j0]

for i in consumers.index:
    # if i != consumers.index[-1]:
    if i != consumers.index[0]:
        # All but one consumer define their mass
        j_f = pp.create_junction(
            net, pn_bar=pn_bar, tfluid_k=293.15, name=f"junction {i} forward",
            geodata=(10*(i+1), 30))
        j_i1 = pp.create_junction(
            net, pn_bar=pn_bar, tfluid_k=293.15, name=f"junction {i} intermediate",
            geodata=(10*(i+1), 20))
        j_i2 = pp.create_junction(
            net, pn_bar=pn_bar, tfluid_k=293.15, name=f"junction {i} intermediate",
            geodata=(10*(i+1), 10))
        j_r = pp.create_junction(
            net, pn_bar=pn_bar, tfluid_k=293.15, name=f"junction {i} return",
            geodata=(10*(i+1), 0))

        # I would assume that I could skip the valve if I define the pressure
        # loss in the heat_exchanger, but the results are very different
        # Probably related to https://github.com/e2nIEE/pandapipes/issues/682
        # maybe fixed with 0.14?
        pp.create_heat_exchanger(
            net, from_junction=j_f, to_junction=j_i1,
            qext_w=consumers.loc[i, 'qext_w'],
            loss_coefficient=pressure_loss_coefficient,
            inner_diameter_mm=200,
            )
        pp.create_flow_control(
            net, from_junction=j_i1, to_junction=j_i2,
            controlled_mdot_kg_per_s=consumers.loc[i, 'mdot_kg_per_s'])
        pp.create_valve(
            net, junction=j_i2, element=j_r, et='ju',
            inner_diameter_mm=200,  # pandapipes 0.14
            # diameter_m=0.2,  # pandapipes 0.13
            # loss_coefficient=pressure_loss_coefficient,
            )

        # heat_consumer seems like an alternative, but has no option for
        # a pressure loss coefficient
        # pp.create_heat_consumer(
        #     net, from_junction=j_f, to_junction=j_r,
        #     qext_w=consumers.loc[i, 'qext_w'],
        #     controlled_mdot_kg_per_s=consumers.loc[i, 'mdot_kg_per_s'],
        #     # deltat_k=dT,
        #     treturn_k=None, name=None, index=None,
        #     in_service=True, type="heat_consumer")

        js_f.append(j_f)
        js_r.append(j_r)
    else:
        # One consumer must not define its mass flow, or else the system
        # becomes over-determined ("Matrix is exactly singular"). Instead
        # it receives "pump mass flow" minus "sum of other consumers",
        # which is fine.
        j_f = pp.create_junction(
            net, pn_bar=pn_bar, tfluid_k=293.15, name=f"junction {i} forward",
            geodata=(10*(i+1), 30))
        j_i = pp.create_junction(
            net, pn_bar=pn_bar, tfluid_k=293.15, name=f"junction {i} intermediate",
            geodata=(10*(i+1), 20))
        j_r = pp.create_junction(
            net, pn_bar=pn_bar, tfluid_k=293.15, name=f"junction {i} return",
            geodata=(10*(i+1), 0))

        pp.create_heat_exchanger(
            net, from_junction=j_f,
            to_junction=j_i,
            # to_junction=j_r,
            qext_w=consumers.loc[i, 'qext_w'],
            loss_coefficient=pressure_loss_coefficient,
            inner_diameter_mm=200,
            )
        pp.create_valve(
            net,
            # junction=j_f,
            junction=j_i,
            element=j_r, et='ju',
            inner_diameter_mm=200,  # pandapipes 0.14
            # diameter_m=0.2,  # pandapipes 0.13
            # loss_coefficient=pressure_loss_coefficient,
            )

        # pp.create_heat_consumer(
        #     net, from_junction=j_f, to_junction=j_r,
        #     qext_w=consumers.loc[i, 'qext_w'],
        #     controlled_mdot_kg_per_s=consumers.loc[i, 'mdot_kg_per_s'],
        #     # deltat_k=dT,
        #     treturn_k=None, name=None, index=None,
        #     in_service=True, type="heat_consumer")

        js_f.append(j_f)
        js_r.append(j_r)


# Create all the pipes
for i in range(len(js_f)):
    if i == 0:
        continue
    pp.create_pipe(net, from_junction=js_f[i-1], to_junction=js_f[i],
                   std_type="ISOPLUS_DRE100_STD", length_km=1, k_mm=0.1, loss_coefficient=0,
                   sections=5, text_k=283, name=f"pipe {i-1}-{i} forward")
    pp.create_pipe(net, from_junction=js_r[i-1], to_junction=js_r[i],
                   std_type="ISOPLUS_DRE100_STD", length_km=1, k_mm=0.1, loss_coefficient=0,
                   sections=5, text_k=283, name=f"pipe {i-1}-{i} return")


pp.plotting.simple_plot(net)

# In pandapipes 0.13, 'bidirectional' and 'sequential' worked
# pp.pipeflow(net, mode='bidirectional', iter=1000)
# pp.pipeflow(net, mode='sequential', iter=1000)
# In pandapipes 0.14, only 'hydraulics' seems to work
pp.pipeflow(net, mode='hydraulics', iter=1000)

df_j = pd.concat([net.junction, net.res_junction, net.junction_geodata], axis='columns')
print(df_j)

# Plot the pressure text at the junction locations
fig, ax = plt.subplots()
ax.scatter(x=df_j['x'], y=df_j['y'], s=50)
for x, y, s in zip(df_j['x'], df_j['y'], df_j['p_bar'].round(4)):
    ax.text(x=x, y=y, s=s)
plt.show()

# The plot shows that there can be problematic cases where the return side
# of a consumer can have a higher pressure than the forward side

# print(net.res_pipe[["mdot_from_kg_per_s"]])
# print(net.res_pipe.T)
# print(net.pipe)

# Plot forward and return pressure distribution
fig, ax = plt.subplots()
ax.plot(range(len(js_f)), net.res_junction.loc[js_f, 'p_bar'], label='Forward')
ax.plot(range(len(js_r)), net.res_junction.loc[js_r, 'p_bar'], label='Return')
ax.set_xlabel('Distance [km]')
ax.set_ylabel('Pressure [bar]')
plt.legend()

# Plot forward and return temperature distribution
fig, ax = plt.subplots()
ax.plot(range(len(js_f)), net.res_junction.loc[js_f, 't_k'], label='Forward')
ax.plot(range(len(js_r)), net.res_junction.loc[js_r, 't_k'], label='Return')
ax.set_xlabel('Distance [km]')
ax.set_ylabel('Temperature [K]')
plt.legend()