A new analytical method for modeling network hydraulics in district systems with bidirectional mass and energy flows

Abstract

Recent technologies of district heating and cooling (DHC) systems integrate renewable energy sources and connect multiple heat producers and consumers to the same network. The resulting interactions may therefore cause the fluid in the pipe, or the energy carried by the system, to flow in both directions. This article presents a newly developed analytical method for evaluating the hydraulic states in DHC networks with bidirectional mass and energy flows. The formulation of the method is presented in detail and is mainly based on the definition of two functions, namely, the flow-pressure function and the bifurcation function. The former describes the relationship between the pressure drop and the mass flow rate in a circular pipe segment, whereas the latter finds the sum of all pressure differences in a closed loop and has a single unique zero that determines the unknown mass flow rate to fulfill network mass balance. Several examples are provided to demonstrate the applicability of the method to single- and multi-looped network topologies. The method enables robust and fast calculations for virtually any network configuration with any number of connected prosumers. It can be easily implemented in any design and simulation tool to account for steady-state or dynamic situations.

Introduction

District heating and cooling (DHC) systems are a key technology for decarbonizing the building sector (H. Lund et al. 2010). They enable the integration of renewable energy sources (RES) to meet European Union (EU) and national targets for improving energy efficiency and reducing carbon emissions (European Commission 2016; Government Offices of Sweden 2016; Government Offices of Sweden 2018; Henrik Lund et al. 2018; Mathiesen et al. 2015; Yuan et al. 2021). The share of different fuel inputs and energy sources for DHC systems varies from one country to another. Sweden, for example, predominantly uses bioenergy and municipal solid waste in the production of district heating (DH) with a reported share of about 70% of the total fuel input in 2020 (Abugabbara, Gehlin, et al. 2023). Denmark, on the other hand, currently has a share of about 20% of RES in the DH system and a framework is already outlined to reach a share of 100% by 2060 (H. Lund et al. 2010). Currently, 43% of the total European DH production comes from renewables (Euroheat & Power 2023). Although DH supplies only 12 and 8% of the European and global final heating in buildings and industry, its production rate increased by about 3% in 2021 and is expected to continue to increase (Euroheat & Power 2023; IEA 2022).

Globally, district cooling (DC) systems deliver much smaller amounts compared to DH, with an estimated amount of 83 TWh (Werner 2017). Two-thirds of the global DC deliveries are geographically located in the Middle East. DC is anticipated to grow in Europe due to hot summer events, with an expected increase of 55–60% in the installed cooling capacity in buildings compared to 2010 (Hitchin, Pout, and Riviere 2013). DC was introduced in Sweden in the 1990s and currently exists in 42 cities with a total network length of 660 km and a delivery of 1 TWh (Energiföretagen Sverige 2022; Energimyndigheten and SCB 2020). The DC system located in Stockholm is one of the largest in Europe with 270 MW of installed capacity in 2020, total delivery of 335 GWh, and a pipe length of 250 km (Jangsten 2020). Cooling is produced by three sources: heat pumps used in the DH system (55%), compression and absorption chillers (27%), and seawater free cooling (18%), with accounted distribution losses of 8% (Stockholm Exergie 2020).

Traditionally, DHC systems consist mainly of three components, namely, (1) a centralized plant, (2) a distribution piping network, and (3) substations to exchange heat with the end users (Frederiksen and Werner 2013). In DH systems, it is common to produce heat with high temperatures in the central plant (Werner 2017). A piping network distributes the heat from the plant to the end customers through a heat-carrying fluid. Substations at the end users transfer the heat into the building through heat exchangers. The heat-carrying fluid in the pipes is restricted to flow in one direction, from the plant toward the end customer, because of the central pumping station. A similar concept is applied to DC systems but for very-low-temperature production (Østergaard et al. 2022). In both systems, the temperature of the medium leaving the plant must meet the requirements of all connected customers. Consequently, distribution losses (or gains) are high and may account for up to one-third of the total delivered energy (Statistics Sweden 2019).

The development of new technologies of DH systems has been centered on lowering the network temperature, which has paved the way toward realizing smart and innovative solutions (Brand et al. 2014; Henrik Lund et al. 2014). As the network operating temperature has been reduced to levels very close to the ambient (e.g., 5–20 °C), the exploitation of low-enthalpy renewable sources like geothermal, solar thermal, and ambient air or water with large shares has become possible. Additionally, multiple low-grade waste heat sources such as sewage networks come into play. The integration of multiple heat sources into the same network poses a technical challenge where the medium direction is no longer restricted by a central plant. Rather, any building connected to the network may become a producer or a consumer of heat (the term prosumer is used in the remainder of the article) (Abugabbara et al. 2020; Brange, Englund, and Lauenburg 2016).

Several network topologies can be adopted for the prosumer concept, including one-, two-, or even three-pipe systems (Arabkoohsar and Alsagri 2020; Buffa et al. 2019; Bünning et al. 2018; Lindhe et al. 2022; Sommer et al. 2020; Zarin Pass, Wetter, and Piette 2018). The most common application found in literature is for two-pipe systems with warm and cold sides. In such a configuration, decentralized heat pumps or chillers extract water from the respective warm and cold pipes to modulate the network temperature to the desired building supply levels. Heat exchangers are also installed at the building level for possible direct supply of cooling and heating depending on the network temperature and the building temperature requirements. Due to the connection of multiple heat sources and depending on whether heating or cooling demands are dominant, the fluid can flow in any direction. However, it is important to distinguish between bidirectional mass and energy flows (Lindhe et al. 2022; Sulzer et al. 2021). The former can be typically observed between seasons when the same building has heating demands in winter and cooling demands in summer and when the fluid extraction changes from the warm to the cold pipe. In the latter case, the type of the dominant demand across the network determines whether the energy flows into or out of the network, regardless of the season. Such bidirectional networks with low temperatures are referred to as bidirectional low temperature networks (BLTNs) (Bu et al. 2023; Licklederer et al. 2021; Wirtz et al. 2020).

The distribution pipes in BLTNs are typically uninsulated due to the low temperature of the heat-carrying fluid. Abugabbara, Javed, and Johansson (2022) presented a Modelica model for steady-state heat loss calculation for uninsulated double buried pipes using the multipole method originally developed by Wallentén (1991) and the equivalent resistance model described by van der Heijde, Aertgeerts, et al. (2017). The Modelica model also accounts for pressure drops due to wall friction that have a quadratic relation to the mass flow rate and includes a constant flow resistance obtained from nominal conditions. The requirement of predefined nominal parameters before solving the network can create a tedious task when modeling BLTNs with multiple prosumers. A thorough hydraulic modeling of BLTNs is therefore needed, given that prosumers extract the heat-carrying fluid at varying amounts and directions, which may also cause circulation pumps to have hydraulic imbalances.

As more prosumers with varying loads are connected to BLTNs, nonlinear hydraulic connections may become apparent. In traditional DHC networks for which the topology is described by a tree, no closed loops are created and hence a linear hydraulic system is formed (Cantoni et al. 2007). On the other hand, BLTNs are initially structured as a one closed-loop network, which after connecting multiple prosumers can become a meshed network consisting of several closed loops (von Rhein et al. 2019). Pressure cycles occur at the connection points between multiple loops and hence a nonlinear hydraulic system is found and the electric pumping power is consequently impacted (De Persis and Kallesoe 2011). It is worth noting that in conventional DHC systems with central plants, pumping energy constitutes less than 0.5% of the total delivered energy (Buffa et al. 2019). However, the share is expected to increase in BLTNs due to the low temperature difference that yields higher flows, as the pumping power is proportional to the third power of the flow.

Much of the research conducted on bidirectional networks investigated the systems using high-level modeling tools incorporating numerical methods (Bilardo et al. 2021; Bu et al. 2023; von Rhein et al. 2019). Abugabbara et al. (2021) compared the pumping energy consumption in two substation design alternatives. The authors reported that adding a free-cooling heat exchanger reduced the total electricity consumption but nearly doubled the share of pumping energy. Wang, Wang, and Zhu (2017) presented a hydraulic regulation method for district systems with distributed variable-speed pumps. Compared to central circulating pumps, the new regulation method was reported to save 26–90% of pumping energy, depending on the rounds of regulation and the frequency adjustment resolution. Sommer et al. (2022) developed a method to alternate the connection of the expansion vessel located at the suction side of the circulation pump. The developed method ensures that enough pressure is maintained to prevent cavitation, and the authors managed to reduce the network pressure by about 8%.

Previous studies on modeling hydraulics of DHC systems and BLTNs utilized numerical methods. Hirsch and Nicolai (2022) presented a thermohydraulic model for a novel DHC network. The hydraulic equations were solved numerically using an embedded Newton–Raphson method and the model was reported to be approximately 75 times faster than a similar model in Modelica. The Modelica language utilizes different integration algorithms to solve the model’s system of equations and has been used in several studies to model district network hydraulics (Abugabbara, Javed, and Johansson 2022; Bünning et al. 2018; Schweiger et al. 2017; Sommer et al. 2020; van der Heijde, Fuchs, et al. 2017). Numerical methods inevitably have common problems that may arise when modeling network hydraulics, such as discretization of time steps, iterations, and convergence errors to approximate the solution. These problems may cause the simulation to halt without providing useful information about the error source.

To support the successful implementation of bidirectional DHC systems, this article presents a new fast analytical method to model network hydraulics for BLTNs with any number of connected prosumers. This work derives from a preliminary study in the report of Lindhe (2022). The method is robust in the sense that it enables modeling different network topologies and hydraulic configurations since a detailed description of the district system is not needed. Additionally, the method provides short modeling and simulation time, which would contribute to an accurate and safe hydraulic design. The newly developed method is implemented using Mathcad software and can be easily deployed to any other software and programming language.

General methods

The article begins in the section “A descriptive network example with three prosumers” with an introduction of the sign convention and notations used in the remainder of the article, with the help of a descriptive BLTN example with three connected prosumers. The section “Flow-pressure formula for any pipe segment” presents the mathematical derivation of the flow-pressure function for different flow regimes where the common equations to describe the pressure drops in a circular pipe segment are rewritten in a unified way. The section “General case with any number of prosumers” introduces the bifurcation function, along with its general expression for any number of connected prosumers, and utilizes its unique zero solution in fulfilling the network mass balance. The section “Turbulent flow in all pipe segments” analyzes the use of the bifurcation function in the turbulent flow regime and provides a simpler equation to solve the unknown mass flow rate in the network. The section “Modifications for two viscosities ${\mu }_{+}=\mu (T_{+})$ and ${\mu }_{-}=\mu (T_{-})$” examines the impact of varying the medium viscosity and reformulates the ordinary equations to determine the node pressure as a function of the viscosity. The section “Two closed loops, two bifurcation functions” investigates the special case of network expansion when a two-closed-loop network is formed. The discussion and conclusions section provides a summary of the full solution, investigates the impact of pipe roughness by implementing another equation to calculate the friction factor, lists the model limitations and briefly describes the comparative validation against three numerical methods, and finally outlines the main conclusions from the study.

A descriptive network example with three prosumers

This section introduces the sign convention and notations used in the remainder of the article with the aid of a network example shown in Figure 1. The network consists of three prosumers that deliver and receive the mass flows ${\dot{M}}_{\nu }$ (kg/s), with $\nu =1,2,3.$ The loop of pipes is connected to an accumulator tank with temperature stratification.

Fig. 1. Piping grid and notations for three prosumers and an accumulator.

All pipe couples have a warm or tepid pipe shown in red, and a cold pipe shown in blue. The cooling or heating demand from each prosumer determines the required mass flows ${\dot{M}}_{\nu }.$ These mass flows may be positive or negative (or zero) to include the bidirectional aspect with positive and negative signs. For positive ${\dot{M}}_{\nu },$ the water flows from the prosumer to the loop in the red pipe, and in the opposite direction in the blue pipe following the arrows in Figure 1. For negative ${\dot{M}}_{\nu },$ these flow directions are reversed.

The closed piping loop consists of $3+1$ pipe segments and nodes (or pipe junctions) enumerated counterclockwise from $\nu =0$ to $\nu =3.$ Node $\nu =0$ connects the accumulator to the loop, and nodes $\nu =1,2,3$ connect the loop to prosumer $\nu .$ Pipe segment $\nu$ lies between node $\nu$ and node $\nu +1$ for $\nu =0,1,$ and $2.$ The last pipe segment $\nu =3$ lies between node 3 and accumulator node $\nu =0.$ Pipe segment $\nu$ has a length $L_{\nu }$ and a hydraulic pipe diameter $D_{\nu }.$

The mass flow process is driven by the prescribed prosumer flows ${\dot{M}}_{\nu }$ with opposite flow directions on the warm (red) and cold (blue) sides. The flow process is exactly antisymmetric between the warm and cold sides. By antisymmetric flow process we mean that the mass flows in the warm and cold sides are equal in magnitude but opposite in direction. However, the fluid viscosity $\mu =\mu (T)$ differs on the two sides due to differences in the fluid temperature. The analysis in this article is first performed for a mean viscosity of warm and cold fluids. The impact of varying viscosities between the warm and cold sides is analyzed in the section “Modifications for two viscosities ${\mu }_{+}=\mu (T_{+})$ and ${\mu }_{-}=\mu (T_{-}).$”

The mass flow in the two pipes of segment $\nu$ in the loop is ${\dot{m}}_{\nu }.$ For positive ${\dot{m}}_{\nu },$ the flow is clockwise in the warm pipe and counterclockwise in the cold pipe, as indicated by the arrows along the pipe loop. For negative ${\dot{m}}_{\nu },$ the flows are reversed.

The mass balances at the nodes connected to the three prosumers and the node connected to the accumulator become (1) $$\begin{array}{c}{\dot{m}}_{\nu }+{\dot{M}}_{\nu }={\dot{m}}_{\nu -1},\hspace{1em}\nu =1,2,3;\hspace{1em}{\dot{m}}_1={\dot{m}}_0-{\dot{M}}_1,\hspace{1em}{\dot{m}}_2={\dot{m}}_1-{\dot{M}}_2={\dot{m}}_0-{\dot{M}}_1-{\dot{M}}_2,\\ \\ {\dot{m}}_3={\dot{m}}_2-{\dot{M}}_3={\dot{m}}_0-{\dot{M}}_1-{\dot{M}}_2-{\dot{M}}_3;\hspace{1em}{\dot{m}}_0-{\dot{m}}_3={\dot{M}}_1+{\dot{M}}_2+{\dot{M}}_3\end{array}$$(1)

The mass flow through the accumulator (positive for downward flow and negative for upward flow) is equal to the sum of the prosumer flows. The mass balances at all nodes are fulfilled for any value of ${\dot{m}}_0.$ The determination of all mass flows ${\dot{m}}_{\nu }$ requires one more relation.

The excess pressure $p_0$ above the accumulator pressure at node $\nu =0$ is zero. The (excess) pressures on the cold side of the nodes then become $-p_{\nu },$ since the process is antisymmetric between warm and cold pipes (for a mean viscosity for warm and cold sides). See Figure 1.

The water pressure $p_{\text{acc}}$ at the accumulator is equal to the sum of atmospheric pressure $p_{\text{atm}}$ and a suitable overpressure $p_{\mathrm{o}\text{.p}.}.$ Let $p_{\nu }$ denote the pressure above $p_{\text{acc}}$ at node $\nu$ on the warm side: (2) $$p_{\text{acc}}=p_{\text{atm}}+p_{\mathrm{o}\text{.p}.},\hspace{1em}{p}_{\text{total,\hspace{0.17em}warm\hspace{0.17em}\hspace{0.17em}side,}\nu }=p_{\text{acc}}+p_{\nu };\hspace{1em}{p}_0=0$$(2)

Flow-pressure formula for any pipe segment

This section presents the mathematical derivation of the flow-pressure formula $f_p(x)$ and describes its properties. Figure 2 shows a pipe with flowing water. The length of the considered pipe is $L$ and the inner diameter is $D.$ The pressure drop over the pipe is $\Delta p,$ and the pressure level at the right-hand outlet is $p_L.$ The bulk water velocity is $u$ (m/s) and the corresponding mass flow is $\dot{m}$ (kg/s). The water density is $\rho$ (kg/m³) and the viscosity is $\mu$ (N·s/m² = Pa·s = kg/(m·s)). The water flow may be positive, zero, or negative to account for flow reversal. The velocity and mass flow become negative for negative $\Delta p.$

Fig. 2. Notations for a pipe segment.

The viscosity depends on the fluid temperature, which is different on the warm and cold sides of the network. The analysis is first made for a mean viscosity. The necessary modifications to account for the two temperature levels are presented in the section “Modifications for two viscosities ${\mu }_{+}=\mu (T_{+})$ and ${\mu }_{-}=\mu (T_{-}).$”

The mass flow is proportional to the bulk velocity. In the formulas here it is convenient to allow the Reynolds number to be negative for negative $u:$ (3) $$\dot{m}=\rho \cdot \pi \cdot {\left(\frac{D}{2}\right)}^2\cdot u\mathrm{(}\text{kg/s}\mathrm{)},\hspace{1em}\mathrm{\text{Re}}=\frac{\rho Du}{\mu },\hspace{1em}{\mathrm{\text{Re}}}_{\mathit{\text{tr}}}=2000$$(3)

The transition from laminar to turbulent flow occurs around $\left|\mathrm{\text{Re}}\right|={\mathrm{\text{Re}}}_{\mathit{\text{tr}}}=2000.$

Flow-pressure formulas for laminar and turbulent flow

Flow-pressure formulas for laminar and turbulent flow in a circular pipe may be found in many books on mass transfer, for example, Kay and Nedderman (1974). The Darcy–Weisbach equation used to calculate the pressure loss in pipes due to friction in the pipe wall is expressed as (4) $$\Delta p=c_f\cdot \frac{L}{D}\cdot \frac{\rho \cdot u^2}{2}$$(4) where $c_f$ is the pipe friction factor. In the laminar region when the magnitude of the Reynolds number is smaller than the transition value of around 2000 (Shashi Menon 2015), the Hagen–Poiseuille equation (Reay, Kew, and McGlen 2014) is used to describe the friction factor in a circular smooth pipe as (5) $$c_f=\frac{64}{\mathrm{\text{Re}}}$$(5)

Substituting Equation 5(5) $$c_f=\frac{64}{\mathrm{\text{Re}}}$$(5) in Equation 4(4) $$\Delta p=c_f\cdot \frac{L}{D}\cdot \frac{\rho \cdot u^2}{2}$$(4) gives (6) $$\left|\mathrm{\text{Re}}\right|<{\mathrm{\text{Re}}}_{\mathit{\text{tr}}}:\hspace{1em}\Delta p=\frac{32\mu Lu}{D^2}=\frac{32\mu L}{D^2}\cdot \underset{=u}{\underbrace{\frac{\mu \mathrm{\text{Re}}}{\rho D\mathrm{ }}}}=\frac{32{\mu }^2}{\rho }\cdot \frac{L}{D^3}\cdot \mathrm{\text{Re}}$$(6)

The Blasius formula (Massey 2006) is often used for the turbulent flow in a circular, smooth pipe, when the magnitude of Reynolds number exceeds the transition value of 2000, and is expressed as (7) $$\left|\mathrm{\text{Re}}\right|>{\mathrm{\text{Re}}}_{\mathit{\text{tr}}}:\hspace{1em}\Delta p=2c_f\cdot \frac{L}{D}\cdot \rho u\left|u\right|,\hspace{1em}{c}_f=\frac{0.079}{{\left|\mathrm{\text{Re}}\right|}^{0.25}}$$(7)

It should be noted that the normally occurring factor $u^2$ is replaced by $u\left|u\right|$ in order to include the case of negative $u$ values.

The preceding formula may be rewritten as follows: (8) $$\Delta p=\frac{2\cdot 0.079}{{\left|\mathrm{\text{Re}}\right|}^{0.25}}\cdot \frac{L}{D}\cdot \rho \cdot \mathrm{ }{\left(\frac{\mu }{\rho D}\right)}^2\cdot \mathrm{\text{Re}}\mathrm{ }\cdot \left|\mathrm{\text{Re}}\right|=\frac{32{\mu }^2}{\rho }\cdot \frac{L}{D^3}\cdot \frac{2\cdot 0.079}{32}\cdot \mathrm{\text{Re}}\cdot {\left|\mathrm{\text{Re}}\right|}^{0.75}$$(8)

A “normalized” Reynolds number $x$ and the constants $k_0$ and $k_1$ are now introduced with their corresponding units: (9) $$x=\frac{\mathrm{\text{Re}}}{{\mathrm{\text{Re}}}_{\mathit{\text{tr}}}}(‐),\hspace{1em}{k}_0=\frac{32{\mu }^2}{\rho }\cdot {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}\text{(N)},\hspace{1em}{k}_1=\frac{0.079}{16}\cdot {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}^{0.75}\mathrm{(}‐\mathrm{)}$$(9)

The flow is laminar for $-1<x<1,$ and it is turbulent for $x>1$ and $x<-1.$ The following combined formula is obtained from Equations 6(6) $$\left|\mathrm{\text{Re}}\right|<{\mathrm{\text{Re}}}_{\mathit{\text{tr}}}:\hspace{1em}\Delta p=\frac{32\mu Lu}{D^2}=\frac{32\mu L}{D^2}\cdot \underset{=u}{\underbrace{\frac{\mu \mathrm{\text{Re}}}{\rho D\mathrm{ }}}}=\frac{32{\mu }^2}{\rho }\cdot \frac{L}{D^3}\cdot \mathrm{\text{Re}}$$(6) , 8(8) $$\Delta p=\frac{2\cdot 0.079}{{\left|\mathrm{\text{Re}}\right|}^{0.25}}\cdot \frac{L}{D}\cdot \rho \cdot \mathrm{ }{\left(\frac{\mu }{\rho D}\right)}^2\cdot \mathrm{\text{Re}}\mathrm{ }\cdot \left|\mathrm{\text{Re}}\right|=\frac{32{\mu }^2}{\rho }\cdot \frac{L}{D^3}\cdot \frac{2\cdot 0.079}{32}\cdot \mathrm{\text{Re}}\cdot {\left|\mathrm{\text{Re}}\right|}^{0.75}$$(8) , and 9(9) $$x=\frac{\mathrm{\text{Re}}}{{\mathrm{\text{Re}}}_{\mathit{\text{tr}}}}(‐),\hspace{1em}{k}_0=\frac{32{\mu }^2}{\rho }\cdot {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}\text{(N)},\hspace{1em}{k}_1=\frac{0.079}{16}\cdot {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}^{0.75}\mathrm{(}‐\mathrm{)}$$(9) : (10) $$\Delta p=k_0\cdot \frac{L}{D^3}\cdot f_p(x),\hspace{1em}{f}_p(x)=\{\begin{array}{cc}x\hspace{1em}& \left|x\right|\le 1\\ k_1\cdot x\cdot {\left|x\right|}^{0.75}& \left|x\right|>1\end{array},\hspace{1em}x=\frac{\mathrm{\text{Re}}}{{\mathrm{\text{Re}}}_{\mathit{\text{tr}}}}=\frac{\rho D\mathrm{ }u}{\mu {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}}$$(10)

Graphs for flow-pressure formula f_p(x)

The dimensionless flow-pressure function $f_p(x)$ is shown in Figures 3 and 4 for $k_1=1.477.$ It is an odd function of $x$ that increases monotonously with $x.$

Fig. 3. Flow-pressure function f_p(x).

Fig. 4. Flow-pressure function f_p(x) up to larger x.

The dashed red curve shows the function with its two branches. It follows $x,$ the blue dotted curve for laminar flow, in $-1\le x\le 1.$ Outside the laminar region, $x>1$ and $x<-1,$ it follows the black curve, where $k_1\cdot x\cdot {(\left|x\right|)}^{0.75}$ is the magnitude of the function $f_p(x)$ in the turbulent flow, as previously introduced in Equations 9(9) $$x=\frac{\mathrm{\text{Re}}}{{\mathrm{\text{Re}}}_{\mathit{\text{tr}}}}(‐),\hspace{1em}{k}_0=\frac{32{\mu }^2}{\rho }\cdot {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}\text{(N)},\hspace{1em}{k}_1=\frac{0.079}{16}\cdot {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}^{0.75}\mathrm{(}‐\mathrm{)}$$(9) and 10(10) $$\Delta p=k_0\cdot \frac{L}{D^3}\cdot f_p(x),\hspace{1em}{f}_p(x)=\{\begin{array}{cc}x\hspace{1em}& \left|x\right|\le 1\\ k_1\cdot x\cdot {\left|x\right|}^{0.75}& \left|x\right|>1\end{array},\hspace{1em}x=\frac{\mathrm{\text{Re}}}{{\mathrm{\text{Re}}}_{\mathit{\text{tr}}}}=\frac{\rho D\mathrm{ }u}{\mu {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}}$$(10) . The flow function is discontinuous at $x=1$ and $x=-1$ with the jump $k_1-1.$ For $x>1$ the function is proportional to $x$ to the power 1.75.

The flow in any pipe segment will typically be turbulent in the applications considered in this study. The difference if the turbulent formula is used in the laminar region becomes (11) $$f_{p,\text{turb}}(x)=k_1\cdot x\cdot {\left|x\right|}^{0.75},\hspace{1em}{f}_p(x)-f_{p,\text{turb}}(x)=\{\begin{array}{cc}x-k_1\cdot x\cdot {\left|x\right|}^{0.75}& \left|x\right|\le 1\\ 0\hspace{1em}& \left|x\right|>1\end{array}$$(11)

Figure 5 shows this difference $f_p(x)-f_{p,\text{turb}}(x).$ The difference varies between 0 and 0.477, which occur respectively at $\left|x\right|=0.5945$ and exactly at transition at $x=1$ and $x=-1.$ Accordingly, the obtained pressure drop can differ by 0% up to 47.7% if the turbulent formula is used in the laminar region.

Fig. 5. The difference $f_p(x)-f_{p\_\text{turb}}(x).$

The flow-pressure formula in Equation 10(10) $$\Delta p=k_0\cdot \frac{L}{D^3}\cdot f_p(x),\hspace{1em}{f}_p(x)=\{\begin{array}{cc}x\hspace{1em}& \left|x\right|\le 1\\ k_1\cdot x\cdot {\left|x\right|}^{0.75}& \left|x\right|>1\end{array},\hspace{1em}x=\frac{\mathrm{\text{Re}}}{{\mathrm{\text{Re}}}_{\mathit{\text{tr}}}}=\frac{\rho D\mathrm{ }u}{\mu {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}}$$(10) gives the pressure drop as a function of $x$ (or $u$). An alternative is to use $x$ as function of the mass flow $\dot{m}:$ (12) $$\dot{m}=\frac{\pi \rho D^2\cdot u}{4}=\frac{\pi \mu D\cdot \mathrm{\text{Re}}}{4},\hspace{1em}x=\frac{\mathrm{\text{Re}}}{{\mathrm{\text{Re}}}_{\mathit{\text{tr}}}}=\frac{4\dot{m}}{\pi \mu D\cdot {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}}=\frac{\dot{m}}{k_2\cdot D};\hspace{1em}{k}_2=\frac{\pi \mu {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}}{4}$$(12)

The new constant $k_2$ has the same dimension as the viscosity: kg/(m·s).

Now three constants are used: (13) $$k_0=\frac{32{\mu }^2}{\rho }\cdot {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}\text{(N)},\hspace{1em}{k}_1=\frac{0.079}{16}\cdot {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}^{0.75}\mathrm{(}‐\mathrm{)},\hspace{1em}{k}_2=\frac{\pi \mu {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}}{4}\mathrm{(}\text{kg/}\mathrm{(}\mathrm{m}\cdot \text{s))}$$(13)

The following input data has been used when generating Figures 3–5: (14) $$\rho =1000(\mathrm{\text{kg}}/{\mathrm{m}}^3),\hspace{1em}\mu =\mu ({20}^{\mathrm{o}}\mathrm{C})=1.00\cdot {10}^{-3}(\mathrm{\text{kg}}/\mathrm{}(\mathrm{m}\cdot \mathrm{s})),\hspace{1em}{\mathrm{\text{Re}}}_{\mathit{\text{tr}}}=2000$$(14)

The three constants of Equation 13(13) $$k_0=\frac{32{\mu }^2}{\rho }\cdot {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}\text{(N)},\hspace{1em}{k}_1=\frac{0.079}{16}\cdot {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}^{0.75}\mathrm{(}‐\mathrm{)},\hspace{1em}{k}_2=\frac{\pi \mu {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}}{4}\mathrm{(}\text{kg/}\mathrm{(}\mathrm{m}\cdot \text{s))}$$(13) become (15) $$k_0=64\cdot {10}^{-6}(\mathrm{N}),\hspace{1em}{k}_1=1.477\mathrm{(}‐\mathrm{)},\hspace{1em}{k}_2=1.571\mathrm{(}\text{kg/}\mathrm{(}\mathrm{m}\cdot \text{s))}$$(15)

General case with any number of prosumers

This section presents the equations to calculate the mass flows ${\dot{m}}_{\nu }$ for any number of prosumers. The key so-called bifurcation function is introduced and discussed in some detail. Figure 6 shows the warm side of the network for any number ${\nu }_{\mathit{\text{ps}}}$ of prosumers. The cold part may be omitted due to the antisymmetry of the flow process. It is worth noting that any prosumer flow ${\dot{M}}_{\nu }$ may in actual cases represent the sum of several branches in the district system, as shown in the magnified part in Figure 6. The mass flows and the (excess) pressures at nodes are the same as on the warm side but with opposite signs (for a mean viscosity). The piping loop consists of ${\nu }_{\mathit{\text{ps}}}+1$ segments and nodes (or pipe junctions) enumerated counterclockwise from $\nu =0$ to ${\nu }_{\mathit{\text{ps}}}.$

Fig. 6. Network for the pipes on the warm side for ${\nu }_{\mathit{\text{ps}}}$ prosumers.

The mass flow process is driven by the prescribed prosumer flows ${\dot{M}}_{\nu }$ with opposite flow directions on the warm (red) and cold sides. The flow process is exactly antisymmetric between the warm and cold sides, except for the fact that the fluid viscosity $\mu =\mu (T)$ differs on the two sides. The analysis is first performed for a mean viscosity of warm and cold fluids. The modifications to account for different warm and cold viscosities are dealt with in the section “Modifications for two viscosities ${\mu }_{+}=\mu (T_{+})$ and ${\mu }_{-}=\mu (T_{-}).$”

The mass flow in pipe segment $\nu$ in the loop is ${\dot{m}}_{\nu }.$ For positive ${\dot{m}}_{\nu },$ the flow is clockwise in warm pipes, as indicated by the arrows along the pipe loop. For negative ${\dot{m}}_{\nu },$ the flow is reversed.

Mass balances

The mass balances at the nodes connected to the ${\nu }_{\mathit{\text{ps}}}$ prosumers become (16) $$\nu =1,\dots {\nu }_{\text{ps}}:\hspace{1em}{\dot{m}}_{\nu }+{\dot{M}}_{\nu }={\dot{m}}_{\nu -1},\hspace{1em}\sum _{j=1}^{\nu }\left({\dot{m}}_j-{\dot{m}}_{j-1}\right)={\dot{m}}_{\nu }-{\dot{m}}_0=-\sum _{j=1}^{\nu }{\dot{M}}_j$$(16)

Let ${\Sigma }_{\nu }$ denote the right-hand sum of accumulated prosumer flows. The following concise formulas for the flows ${\dot{m}}_{\nu }$ are obtained: (17) $$\boxed{{\dot{m}}_{\nu }={\dot{m}}_0-{\Sigma }_{\nu },\hspace{1em}\nu =0,1,\dots {\nu }_{\mathit{\text{ps}}};\hspace{1em}{\Sigma }_0=0,\hspace{1em}\nu \ge 1:{\Sigma }_{\nu }=\sum _{j=1}^{\nu }{\dot{M}}_j}$$(17)

This formula follows directly from the mass balance for the rectangular region within the dashed closed line in Figure 6. The mass balance for node 0 connected to the accumulator becomes (18) $${\dot{m}}_0-{\dot{m}}_{{\nu }_{\text{ps}}}={\Sigma }_{{\nu }_{\text{ps}}}={\dot{M}}_1+\dots +{\dot{M}}_{{\nu }_{\text{ps}}}$$(18)

The pressure difference over pipe segment $\nu$ with the mass flow ${\dot{m}}_{\nu }$ from Equations 10(10) $$\Delta p=k_0\cdot \frac{L}{D^3}\cdot f_p(x),\hspace{1em}{f}_p(x)=\{\begin{array}{cc}x\hspace{1em}& \left|x\right|\le 1\\ k_1\cdot x\cdot {\left|x\right|}^{0.75}& \left|x\right|>1\end{array},\hspace{1em}x=\frac{\mathrm{\text{Re}}}{{\mathrm{\text{Re}}}_{\mathit{\text{tr}}}}=\frac{\rho D\mathrm{ }u}{\mu {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}}$$(10) , 12(12) $$\dot{m}=\frac{\pi \rho D^2\cdot u}{4}=\frac{\pi \mu D\cdot \mathrm{\text{Re}}}{4},\hspace{1em}x=\frac{\mathrm{\text{Re}}}{{\mathrm{\text{Re}}}_{\mathit{\text{tr}}}}=\frac{4\dot{m}}{\pi \mu D\cdot {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}}=\frac{\dot{m}}{k_2\cdot D};\hspace{1em}{k}_2=\frac{\pi \mu {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}}{4}$$(12) , and 17(17) $$\boxed{{\dot{m}}_{\nu }={\dot{m}}_0-{\Sigma }_{\nu },\hspace{1em}\nu =0,1,\dots {\nu }_{\mathit{\text{ps}}};\hspace{1em}{\Sigma }_0=0,\hspace{1em}\nu \ge 1:{\Sigma }_{\nu }=\sum _{j=1}^{\nu }{\dot{M}}_j}$$(17) becomes (19) $$p_{\nu +1}-p_{\nu }=k_0\cdot \frac{L_{\nu }}{D_{\nu }^3}\cdot f_p\left(\frac{{\dot{m}}_0-{\Sigma }_{\nu }}{k_2{D}_{\nu }}\right),\hspace{1em}\nu =0,1,\dots {\nu }_{\text{ps}};\hspace{1em}{p}_{{\nu }_{\text{ps}}+1}=p_0$$(19)

The pressure difference over the last pipe segment ${\nu }_{\text{ps}}$ is equal to the pressure $p_0$ at the accumulator node minus the pressure at node ${\nu }_{\text{ps}}.$ The formula is applicable to the last pipe segment ${\nu }_{\text{ps}},$ if $p_{{\nu }_{\text{ps}}+1}$ is equal to $p_0.$

The sum of the preceding pressure differences around the closed loop shown in Figure 6 defines a function of the unknown mass flow $\dot{m}={\dot{m}}_0$ in pipe segment $\nu =0:$ (20) $$\sum _{\nu =0}^{{\nu }_{\text{ps}}}\left(p_{\nu +1}-p_{\nu }\right)=\Delta p_{\text{bf}}(\dot{m}),\hspace{1em}\boxed{\Delta p_{\text{bf}}(\dot{m})=k_0\cdot \sum _{\nu =0}^{{\nu }_{\text{ps}}}\frac{L_{\nu }}{D_{\nu }^3}\cdot f_p\left(\frac{\dot{m}-{\Sigma }_{\nu }}{k_2{D}_{\nu }}\right)}$$(20)

This sum of these pressure differences over the closed loop shown in Figure 6 is called the bifurcation function.

The sum of pressure differences over the closed loop shown in Figure 6 becomes zero since the sequence of node pressure starts and ends with the node pressure $p_0$ at the accumulator node: (21) $$\sum _{\nu =0}^{{\nu }_{\text{ps}}}\left(p_{\nu +1}-p_{\nu }\right)=p_{{\nu }_{\text{ps}}+1}-p_0=0\hspace{1em}\Rightarrow \hspace{1em}\boxed{\Delta p_{\text{bf}}({\dot{m}}_0)=0}$$(21)

The bifurcation function $\Delta p_{\text{bf}}(\dot{m})$ is equal to zero for $\dot{m}={\dot{m}}_0.$ As previously noted, the flow-pressure function $f_p(x)$ is an odd function that increases monotonously for increasing $x.$ This means that the sum $\Delta p_{\text{bf}}(\dot{m})$ also increases monotonously with increasing $\dot{m}$ from large negative values and ending with large positive values. The bifurcation function $\Delta p_{\text{bf}}(\dot{m})$ has a single unique zero: $\Delta p_{\text{bf}}({\dot{m}}_0)=0.$

Bifurcation function for a single prosumer

Figure 7 shows the network for the simple case of a single prosumer with a (positive) prosumer flow ${\dot{M}}_1.$ The part ${\dot{m}}_0$ of ${\dot{M}}_1$ flows along the upper pipe (pipe segment $\nu =0$), and the remaining part ${\dot{M}}_1-{\dot{m}}_0=-{\dot{m}}_1$ flows along the lower pipe (segment $\nu =1$). There is a bifurcation flow at the prosumer node $\nu =1.$ The two flows ${\dot{m}}_0$ and ${\dot{M}}_1-{\dot{m}}_0$ meet at the accumulator node $\nu =0$ in a confluent outflow to the accumulator.

Fig. 7. Network for a single prosumer. Bifurcation flow from the prosumer (right).

The right-hand sum in Equation 20(20) $$\sum _{\nu =0}^{{\nu }_{\text{ps}}}\left(p_{\nu +1}-p_{\nu }\right)=\Delta p_{\text{bf}}(\dot{m}),\hspace{1em}\boxed{\Delta p_{\text{bf}}(\dot{m})=k_0\cdot \sum _{\nu =0}^{{\nu }_{\text{ps}}}\frac{L_{\nu }}{D_{\nu }^3}\cdot f_p\left(\frac{\dot{m}-{\Sigma }_{\nu }}{k_2{D}_{\nu }}\right)}$$(20) defines the bifurcation function, $\Delta p_{\text{bf}}(\dot{m}),$ of the unknown mass flow $\dot{m}:$ (22) $${\nu }_{\text{ps}}=1:\hspace{1em}\Delta p_{\text{bf}}(\dot{m})=k_0\cdot \frac{L_0}{D_0^3}\cdot f_p\left(\frac{\dot{m}-0}{k_2{D}_0}\right)+k_0\cdot \frac{L_1}{D_1^3}\cdot f_p\left(\frac{\dot{m}-{\dot{M}}_1}{k_2{D}_1}\right),\hspace{1em}\Delta p_{\text{bf}}({\dot{m}}_0)=0$$(22)

The root, where the bifurcation function $\Delta p_{\text{bf}}(\dot{m})$ becomes zero, gives the unknown mass flow ${\dot{m}}_0.$

Figure 8 shows this bifurcation function for the case of (23) $$\begin{array}{l}{\nu }_{\text{ps}}=1,\hspace{1em}{\dot{M}}_1=10\mathrm{(}\text{kg/s}\mathrm{)},\hspace{1em}{L}_0=100\text{(m)},\hspace{1em}{L}_1=150\text{(m)},\\ \\ D_0=D_1=0.25\text{(m)},\hspace{1em}\mu =1.0\cdot {10}^{-3}\mathrm{(}\text{kg/}\mathrm{(}\mathrm{m}\cdot \text{s))};\Rightarrow \\ \\ {\Sigma }_0=0,\hspace{1em}{\Sigma }_1={\dot{M}}_1=10;\Delta p_{\text{bf}}(0)<0,\Delta p_{\text{bf}}({\dot{M}}_1)>0\Rightarrow 0<{\dot{m}}_0<10\end{array}$$(23)

Fig. 8. Bifurcation function from Equations 22(22) $${\nu }_{\text{ps}}=1:\hspace{1em}\Delta p_{\text{bf}}(\dot{m})=k_0\cdot \frac{L_0}{D_0^3}\cdot f_p\left(\frac{\dot{m}-0}{k_2{D}_0}\right)+k_0\cdot \frac{L_1}{D_1^3}\cdot f_p\left(\frac{\dot{m}-{\dot{M}}_1}{k_2{D}_1}\right),\hspace{1em}\Delta p_{\text{bf}}({\dot{m}}_0)=0$$(22) and 23(23) $$\begin{array}{l}{\nu }_{\text{ps}}=1,\hspace{1em}{\dot{M}}_1=10\mathrm{(}\text{kg/s}\mathrm{)},\hspace{1em}{L}_0=100\text{(m)},\hspace{1em}{L}_1=150\text{(m)},\\ \\ D_0=D_1=0.25\text{(m)},\hspace{1em}\mu =1.0\cdot {10}^{-3}\mathrm{(}\text{kg/}\mathrm{(}\mathrm{m}\cdot \text{s))};\Rightarrow \\ \\ {\Sigma }_0=0,\hspace{1em}{\Sigma }_1={\dot{M}}_1=10;\Delta p_{\text{bf}}(0)<0,\Delta p_{\text{bf}}({\dot{M}}_1)>0\Rightarrow 0<{\dot{m}}_0<10\end{array}$$(23) with its two parts: $\Delta p_{\text{bf}}(\dot{m})=\Delta p(\dot{m},0)+\Delta p(\dot{m},1).$

The bifurcation function becomes zero at $\dot{m}={\dot{m}}_0=5.577.$ The two terms of the bifurcation function in Equation 22(22) $${\nu }_{\text{ps}}=1:\hspace{1em}\Delta p_{\text{bf}}(\dot{m})=k_0\cdot \frac{L_0}{D_0^3}\cdot f_p\left(\frac{\dot{m}-0}{k_2{D}_0}\right)+k_0\cdot \frac{L_1}{D_1^3}\cdot f_p\left(\frac{\dot{m}-{\dot{M}}_1}{k_2{D}_1}\right),\hspace{1em}\Delta p_{\text{bf}}({\dot{m}}_0)=0$$(22) , $\Delta p(\dot{m},0)$ and $\Delta p(\dot{m},1),$ are defined in Equation 24(24) $$\Delta p_{\text{bf}}(\dot{m})=\sum _{\nu =0}^{{\nu }_{\text{ps}}}\Delta p(\dot{m},\nu ),\hspace{1em}\Delta p(\dot{m},\nu )=k_0\cdot \frac{L_{\nu }}{D_{\nu }^3}\cdot f_p\left(\frac{\dot{m}-{\Sigma }_{\nu }}{k_2{D}_{\nu }}\right);\hspace{1em}{\Sigma }_0=0,\hspace{1em}\nu \ge 1:{\Sigma }_{\nu }=\sum _{j=1}^{\nu }{\dot{M}}_j$$(24) . These two functions balance each other at $\dot{m}={\dot{m}}_0=5.577:$ $-\Delta p({\dot{m}}_0,1)=\Delta p({\dot{m}}_0,0)=62.8.$

An exact formula for solving ${\dot{m}}_0$ for turbulent flows in the case ${\nu }_{\text{ps}}=1$ is given later by Equation 33(33) $${\left(\frac{L_0}{L_1}\cdot \frac{D_1^{4.75}}{D_0^{4.75}}\right)}^{\frac{4}{7}}\cdot {\dot{m}}_0={\dot{M}}_1-{\dot{m}}_0,\hspace{1em}4.75=\frac{19}{4}\hspace{1em}\Rightarrow \hspace{1em}{\dot{m}}_0=\frac{{\dot{M}}_1}{1+{\left(L_0/L_1\right)}^{\frac{4}{7}}\cdot {\left(D_1/D_0\right)}^{\frac{19}{7}}}$$(33) .

Bifurcation function for any number of prosumers

For any number of ${\nu }_{\text{ps}}$ prosumers, the bifurcation function is, from Equation 20(20) $$\sum _{\nu =0}^{{\nu }_{\text{ps}}}\left(p_{\nu +1}-p_{\nu }\right)=\Delta p_{\text{bf}}(\dot{m}),\hspace{1em}\boxed{\Delta p_{\text{bf}}(\dot{m})=k_0\cdot \sum _{\nu =0}^{{\nu }_{\text{ps}}}\frac{L_{\nu }}{D_{\nu }^3}\cdot f_p\left(\frac{\dot{m}-{\Sigma }_{\nu }}{k_2{D}_{\nu }}\right)}$$(20) , (24) $$\Delta p_{\text{bf}}(\dot{m})=\sum _{\nu =0}^{{\nu }_{\text{ps}}}\Delta p(\dot{m},\nu ),\hspace{1em}\Delta p(\dot{m},\nu )=k_0\cdot \frac{L_{\nu }}{D_{\nu }^3}\cdot f_p\left(\frac{\dot{m}-{\Sigma }_{\nu }}{k_2{D}_{\nu }}\right);\hspace{1em}{\Sigma }_0=0,\hspace{1em}\nu \ge 1:{\Sigma }_{\nu }=\sum _{j=1}^{\nu }{\dot{M}}_j$$(24)

Figure 9 shows this bifurcation function for the case of (25) $$\begin{array}{l}{\nu }_{\text{ps}}=5;\hspace{1em}\nu =0,1\dots {\nu }_{\text{ps}}:\hspace{1em}{L}_{\nu }=100+10\cdot \nu \text{(m)},D_{\nu }=0.25\text{(m)};\hspace{1em}\\ \\ \mu =1.0\cdot {10}^{-3}\mathrm{(}\text{kg/}\mathrm{(}\mathrm{m}\cdot \text{s))},\hspace{1em}{k}_2=1.571\mathrm{(}\text{kg/}\mathrm{(}\mathrm{m}\cdot \text{s))};\\ \\ {\dot{M}}_1=10,\hspace{1em}{\dot{M}}_2=-15,\hspace{1em}{\dot{M}}_3=8,\hspace{1em}{\dot{M}}_4=0,\hspace{1em}{\dot{M}}_5=-5\mathrm{(}\text{kg/s}\mathrm{)}\hspace{1em}\Rightarrow \\ \\ {\Sigma }_0=0,\hspace{1em}{\Sigma }_1=10,\hspace{1em}{\Sigma }_2=-5,\hspace{1em}{\Sigma }_3=3,\hspace{1em}{\Sigma }_4=3,\hspace{1em}{\Sigma }_5=-2\end{array}$$(25)

Fig. 9. Bifurcation function $\Delta p_{\text{bf}}(\dot{m}),$ from Equations 24(24) $$\Delta p_{\text{bf}}(\dot{m})=\sum _{\nu =0}^{{\nu }_{\text{ps}}}\Delta p(\dot{m},\nu ),\hspace{1em}\Delta p(\dot{m},\nu )=k_0\cdot \frac{L_{\nu }}{D_{\nu }^3}\cdot f_p\left(\frac{\dot{m}-{\Sigma }_{\nu }}{k_2{D}_{\nu }}\right);\hspace{1em}{\Sigma }_0=0,\hspace{1em}\nu \ge 1:{\Sigma }_{\nu }=\sum _{j=1}^{\nu }{\dot{M}}_j$$(24) and 25(25) $$\begin{array}{l}{\nu }_{\text{ps}}=5;\hspace{1em}\nu =0,1\dots {\nu }_{\text{ps}}:\hspace{1em}{L}_{\nu }=100+10\cdot \nu \text{(m)},D_{\nu }=0.25\text{(m)};\hspace{1em}\\ \\ \mu =1.0\cdot {10}^{-3}\mathrm{(}\text{kg/}\mathrm{(}\mathrm{m}\cdot \text{s))},\hspace{1em}{k}_2=1.571\mathrm{(}\text{kg/}\mathrm{(}\mathrm{m}\cdot \text{s))};\\ \\ {\dot{M}}_1=10,\hspace{1em}{\dot{M}}_2=-15,\hspace{1em}{\dot{M}}_3=8,\hspace{1em}{\dot{M}}_4=0,\hspace{1em}{\dot{M}}_5=-5\mathrm{(}\text{kg/s}\mathrm{)}\hspace{1em}\Rightarrow \\ \\ {\Sigma }_0=0,\hspace{1em}{\Sigma }_1=10,\hspace{1em}{\Sigma }_2=-5,\hspace{1em}{\Sigma }_3=3,\hspace{1em}{\Sigma }_4=3,\hspace{1em}{\Sigma }_5=-2\end{array}$$(25) , with its six parts $\Delta p(\dot{m},\nu ),$ $\nu =0,1,\dots 5.$

The function $\Delta p(\dot{m},\nu ),$ Equation 24(24) $$\Delta p_{\text{bf}}(\dot{m})=\sum _{\nu =0}^{{\nu }_{\text{ps}}}\Delta p(\dot{m},\nu ),\hspace{1em}\Delta p(\dot{m},\nu )=k_0\cdot \frac{L_{\nu }}{D_{\nu }^3}\cdot f_p\left(\frac{\dot{m}-{\Sigma }_{\nu }}{k_2{D}_{\nu }}\right);\hspace{1em}{\Sigma }_0=0,\hspace{1em}\nu \ge 1:{\Sigma }_{\nu }=\sum _{j=1}^{\nu }{\dot{M}}_j$$(24) center, is zero at $\dot{m}={\Sigma }_{\nu }.$ It is an odd function of $\dot{m}-{\Sigma }_{\nu }.$ The order of the six functions $\Delta (\dot{m},\nu )$ along the $\dot{m}$ axes in Figure 9 depends on the values ${\Sigma }_{\nu }:$ (26) $${\Sigma }_2=-5,\hspace{1em}{\Sigma }_5=-2,\hspace{1em}{\Sigma }_0=0,\hspace{1em}{\Sigma }_3=3,\hspace{1em}{\Sigma }_4=3,\hspace{1em}{\Sigma }_1=10$$(26)

All six functions are positive for $\dot{m}>{\Sigma }_1=10$ and negative for $\dot{m}<{\Sigma }_2=-5.$ The root ${\dot{m}}_0$ must lie in the interval $-5<\dot{m}<10.$ Calculation gives the root ${\dot{m}}_0=1.643.$

For the case presented in this section, the excess pressure at each node on the warm side, $p_{\nu },$ with a single mean viscosity $\mu =\mu (20)$ is shown in Figure 10.

Fig. 10. Excess pressure $p_{\nu }$ for the network example with 5 prosumers. The start and end nodes $\nu =0$ and $\nu =6$ represent the accumulator node.

Turbulent flow in all pipe segments

This section analyzes the bifurcation function for turbulent flows in the pipes. The mass flow function $f_p(x)$ for turbulent flow, $\left|x\right|>1,$ reads, from Equation 10(10) $$\Delta p=k_0\cdot \frac{L}{D^3}\cdot f_p(x),\hspace{1em}{f}_p(x)=\{\begin{array}{cc}x\hspace{1em}& \left|x\right|\le 1\\ k_1\cdot x\cdot {\left|x\right|}^{0.75}& \left|x\right|>1\end{array},\hspace{1em}x=\frac{\mathrm{\text{Re}}}{{\mathrm{\text{Re}}}_{\mathit{\text{tr}}}}=\frac{\rho D\mathrm{ }u}{\mu {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}}$$(10) , (27) $$f_{\text{p\_turb}}(x)=k_1\cdot x{\left|x\right|}^{0.75},\hspace{1em}x=\frac{\dot{m}}{k_2\cdot D}\Rightarrow f_{\text{p\_turb}}(\dot{m})=k_1\cdot \frac{\dot{m}\cdot |\dot{m}{|}^{0.75}}{{(k_2\cdot D)}^{1.75}},\hspace{1em}\left|\dot{m}\right|\ge k_2\cdot D$$(27)

For turbulent flow in all pipe segments, the bifurcation function becomes (28) $$\Delta p_{\text{bf\_turb}}(\dot{m})=\sum _{\nu =0}^{{\nu }_{\text{ps}}}{k}_0\cdot \frac{L_{\nu }}{D_{\nu }^3}\cdot k_1\frac{\left(\dot{m}-{\Sigma }_{\nu }\right)\cdot {\left(|\dot{m}-{\Sigma }_{\nu }|\right)}^{0.75}}{{\left(k_2\cdot D_{\nu }\right)}^{1.75}},\hspace{1em}|\dot{m}-{\Sigma }_{\nu }|>k_2\cdot D_{\nu }$$(28)

The dependence of the various parameters is seen more clearly in the formula (29) $$\Delta p_{\text{bf\_turb}}(\dot{m})=\frac{k_0\cdot k_1}{k_2{}^{1.75}}\cdot \sum _{\nu =0}^{{\nu }_{\text{ps}}}\frac{L_{\nu }}{D_{\nu }^{4.75}}\cdot \left(\dot{m}-{\Sigma }_{\nu }\right)\cdot {\left(|\dot{m}-{\Sigma }_{\nu }|\right)}^{0.75}$$(29)

The relation ${\Delta }_{\text{bf\_turb}}({\dot{m}}_0)=0$ gives the following simpler equations to determine ${\dot{m}}_0:$ (30) $$\boxed{\sum _{\nu =0}^{{\nu }_{\text{ps}}}\frac{L_{\nu }}{D_{\nu }^{4.75}}\cdot \left({\dot{m}}_0-{\Sigma }_{\nu }\right)\cdot {\left(|{\dot{m}}_0-{\Sigma }_{\nu }|\right)}^{0.75}=0,{\Sigma }_0=0,\nu \ge 1:{\Sigma }_{\nu }=\sum _{j=1}^{\nu }{\dot{M}}_j}$$(30)

The mass flow ${\dot{m}}_0$ depends on the prosumer flows ${\dot{M}}_{\nu }$ and the factors $L_{\nu }\cdot {(D_{\nu })}^{-4.75}.$

For applications considered in this study, the (nonzero) prosumer flows correspond to turbulent flow: $\left|{\dot{M}}_{\nu }\right|\ge k_2\cdot D_{\nu }$ for all $\nu .$ This means that all or nearly all of the terms in Equation 30(30) $$\boxed{\sum _{\nu =0}^{{\nu }_{\text{ps}}}\frac{L_{\nu }}{D_{\nu }^{4.75}}\cdot \left({\dot{m}}_0-{\Sigma }_{\nu }\right)\cdot {\left(|{\dot{m}}_0-{\Sigma }_{\nu }|\right)}^{0.75}=0,{\Sigma }_0=0,\nu \ge 1:{\Sigma }_{\nu }=\sum _{j=1}^{\nu }{\dot{M}}_j}$$(30) are the exact turbulent ones. A typical example is shown in Figure 9. Here, all six terms lie by a comfortable margin in the turbulent region: (31) $$\mu =0.001\mathrm{(}\text{kg/}\mathrm{(}\mathrm{m}\cdot \text{s))},\hspace{1em}{k}_2=1.571\mathrm{(}\text{kg/}\mathrm{(}\mathrm{m}\cdot \text{s))},\hspace{1em}{D}_{\nu }=0.25\text{(m)}:\hspace{1em}{k}_2\cdot D_{\nu }=0.393$$(31)

The values of the six functions in their respective laminar intervals lie close to zero. The deviations from zero are hardly visible. This difference between laminar and turbulent expressions for the pipe flow in the laminar interval is given by Equation 11(11) $$f_{p,\text{turb}}(x)=k_1\cdot x\cdot {\left|x\right|}^{0.75},\hspace{1em}{f}_p(x)-f_{p,\text{turb}}(x)=\{\begin{array}{cc}x-k_1\cdot x\cdot {\left|x\right|}^{0.75}& \left|x\right|\le 1\\ 0\hspace{1em}& \left|x\right|>1\end{array}$$(11) and is shown in Figure 5.

Equation 30(30) $$\boxed{\sum _{\nu =0}^{{\nu }_{\text{ps}}}\frac{L_{\nu }}{D_{\nu }^{4.75}}\cdot \left({\dot{m}}_0-{\Sigma }_{\nu }\right)\cdot {\left(|{\dot{m}}_0-{\Sigma }_{\nu }|\right)}^{0.75}=0,{\Sigma }_0=0,\nu \ge 1:{\Sigma }_{\nu }=\sum _{j=1}^{\nu }{\dot{M}}_j}$$(30) for turbulent flows may be used in the present applications provided that the (nonzero) prosumers flows lie above the transition level $\left|{\dot{M}}_{\nu }\right|>k_2\cdot D_{\nu }.$

Equation 30(30) $$\boxed{\sum _{\nu =0}^{{\nu }_{\text{ps}}}\frac{L_{\nu }}{D_{\nu }^{4.75}}\cdot \left({\dot{m}}_0-{\Sigma }_{\nu }\right)\cdot {\left(|{\dot{m}}_0-{\Sigma }_{\nu }|\right)}^{0.75}=0,{\Sigma }_0=0,\nu \ge 1:{\Sigma }_{\nu }=\sum _{j=1}^{\nu }{\dot{M}}_j}$$(30) for ${\dot{m}}_0$ may be solved exactly for the simple case of a single prosumer: (32) $${\nu }_{\text{ps}}=1:\hspace{1em}{\dot{M}}_1>0,\hspace{1em}0<{\dot{m}}_0<{\dot{M}}_1\hspace{1em}\Rightarrow \hspace{1em}\frac{L_0}{D_0^{4.75}}\cdot {\left({\dot{m}}_0\right)}^{1.75}-\frac{L_1}{D_1^{4.75}}\cdot {\left({\dot{M}}_1-{\dot{m}}_0\right)}^{1.75}=0$$(32) (33) $${\left(\frac{L_0}{L_1}\cdot \frac{D_1^{4.75}}{D_0^{4.75}}\right)}^{\frac{4}{7}}\cdot {\dot{m}}_0={\dot{M}}_1-{\dot{m}}_0,\hspace{1em}4.75=\frac{19}{4}\hspace{1em}\Rightarrow \hspace{1em}{\dot{m}}_0=\frac{{\dot{M}}_1}{1+{\left(L_0/L_1\right)}^{\frac{4}{7}}\cdot {\left(D_1/D_0\right)}^{\frac{19}{7}}}$$(33)

The formula is also valid for negative ${\dot{M}}_1.$

Modifications for two viscosities ${\mu }_{+}=\mu (T_{+})$ and ${\mu }_{-}=\mu (T_{-})$

This section examines the impact of varying viscosities on the calculated node pressures. The water viscosity depends on the temperature: $\mu =\mu (T).$ See Table 1.

Table 1. Viscosity of water $\mu (T).$

T (^oC)	0	5	10	15	20	25	30	35
$1000\cdot \mu (T)$	1.788	1.518	1.306	1.137	1.0016	0.890	0.797	0.719
$\sqrt[4]{\mu (T)/\mu (20)}$	1.156	1.110	1.069	1.032	1	0.971	0.945	0.920

The pressures $p_{\nu }$ at the nodes are, for turbulent flows, determined from Equations 19(19) $$p_{\nu +1}-p_{\nu }=k_0\cdot \frac{L_{\nu }}{D_{\nu }^3}\cdot f_p\left(\frac{{\dot{m}}_0-{\Sigma }_{\nu }}{k_2{D}_{\nu }}\right),\hspace{1em}\nu =0,1,\dots {\nu }_{\text{ps}};\hspace{1em}{p}_{{\nu }_{\text{ps}}+1}=p_0$$(19) and 27(27) $$f_{\text{p\_turb}}(x)=k_1\cdot x{\left|x\right|}^{0.75},\hspace{1em}x=\frac{\dot{m}}{k_2\cdot D}\Rightarrow f_{\text{p\_turb}}(\dot{m})=k_1\cdot \frac{\dot{m}\cdot |\dot{m}{|}^{0.75}}{{(k_2\cdot D)}^{1.75}},\hspace{1em}\left|\dot{m}\right|\ge k_2\cdot D$$(27) : (34) $$p_{\nu +1}=p_{\nu }+k_0\cdot \frac{L_{\nu }}{D_{\nu }^3}\cdot k_1\cdot \frac{{\dot{m}}_0-{\Sigma }_{\nu }}{k_2{D}_{\nu }}\cdot \frac{{\left|{\dot{m}}_0-{\Sigma }_{\nu }\right|}^{0.75}}{{\left(k_2{D}_{\nu }\right)}^{0.75}},\hspace{1em}\nu =0,1,\dots {\nu }_{\text{ps}};\hspace{1em}{p}_0=0$$(34)

The factor $k_0$ is proportional to ${\mu }^2$ and $k_2$ is proportional to $\mu :$ (35) $$\frac{k_0\cdot k_1}{{\left(k_2\right)}^{1.75}}=\frac{32{\mu }^2}{\rho }\cdot {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}\cdot \frac{0.079}{16}\cdot {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}^{0.75}\cdot {\left(\frac{4}{\pi \mu {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}}\right)}^{1.75}=k_{\mu }\cdot \sqrt[4]{\mu },\hspace{1em}{k}_{\mu }=\frac{2\cdot 0.079}{\rho }\cdot {\left(\frac{4}{\pi }\right)}^{1.75}=\frac{0.241}{\rho }$$(35)

The pressures $p_{\nu }$ at the nodes for a single $\mu (=\mu (20))$ are given by the preceding relations.

The preceding set of recursive relations may be written: (36) $${\tilde{p}}_{\nu }=\frac{p_{\nu }}{\sqrt[4]{\mu }}:\hspace{1em}{\tilde{p}}_{\nu +1}-{\tilde{p}}_{\nu }=k_{.\mu }\cdot \frac{L_{\nu }}{D_{\nu }^{4.75}}\cdot \left({\dot{m}}_0-{\Sigma }_{\nu }\right)\cdot {|{\dot{m}}_0-{\Sigma }_{\nu }|}^{0.75},\hspace{1em}\nu =0,1,\dots {\nu }_{\text{ps}};\hspace{1em}{\tilde{p}}_0=0$$(36)

The differences ${\tilde{p}}_{\nu +1}-{\tilde{p}}_{\nu }$ and the ensuing sequence of pressure factors ${\tilde{p}}_0,{\tilde{p}}_1,{\tilde{p}}_2\dots$ from Equation 36(36) $${\tilde{p}}_{\nu }=\frac{p_{\nu }}{\sqrt[4]{\mu }}:\hspace{1em}{\tilde{p}}_{\nu +1}-{\tilde{p}}_{\nu }=k_{.\mu }\cdot \frac{L_{\nu }}{D_{\nu }^{4.75}}\cdot \left({\dot{m}}_0-{\Sigma }_{\nu }\right)\cdot {|{\dot{m}}_0-{\Sigma }_{\nu }|}^{0.75},\hspace{1em}\nu =0,1,\dots {\nu }_{\text{ps}};\hspace{1em}{\tilde{p}}_0=0$$(36) are independent of the viscosity. Equation 36(36) $${\tilde{p}}_{\nu }=\frac{p_{\nu }}{\sqrt[4]{\mu }}:\hspace{1em}{\tilde{p}}_{\nu +1}-{\tilde{p}}_{\nu }=k_{.\mu }\cdot \frac{L_{\nu }}{D_{\nu }^{4.75}}\cdot \left({\dot{m}}_0-{\Sigma }_{\nu }\right)\cdot {|{\dot{m}}_0-{\Sigma }_{\nu }|}^{0.75},\hspace{1em}\nu =0,1,\dots {\nu }_{\text{ps}};\hspace{1em}{\tilde{p}}_0=0$$(36) may be applied for a single viscosity as well as for the viscosity on the warm and cold sides: (37) $$p_{\nu }=\sqrt[4]{\mu }\cdot {\tilde{p}}_{\nu },\hspace{1em}{p}_{\nu }^{+}=\sqrt[4]{\mu (T_{+})}\cdot {\tilde{p}}_{\nu },\hspace{1em}{p}_{\nu }^{-}=-\sqrt[4]{\mu (T_{-})}\cdot {\tilde{p}}_{\nu },\hspace{1em}\nu =0,1,\dots {\nu }_{\text{ps}}$$(37)

These three pressure sequences are proportional to the fourth root of the viscosity in accordance with Equation 36(36) $${\tilde{p}}_{\nu }=\frac{p_{\nu }}{\sqrt[4]{\mu }}:\hspace{1em}{\tilde{p}}_{\nu +1}-{\tilde{p}}_{\nu }=k_{.\mu }\cdot \frac{L_{\nu }}{D_{\nu }^{4.75}}\cdot \left({\dot{m}}_0-{\Sigma }_{\nu }\right)\cdot {|{\dot{m}}_0-{\Sigma }_{\nu }|}^{0.75},\hspace{1em}\nu =0,1,\dots {\nu }_{\text{ps}};\hspace{1em}{\tilde{p}}_0=0$$(36) . The three sequences are related to each other by the factor on the third line in Table 1: (38) $$p_{\nu }^{+}=\sqrt[4]{\frac{\mu (T_{+})}{\mu }}\cdot p_{\nu },\hspace{1em}{p}_{\nu }^{-}=-\sqrt[4]{\frac{\mu (T_{-})}{\mu }}\cdot p_{\nu },\hspace{1em}\nu =0,1,\dots {\nu }_{\text{ps}}$$(38)

The temperatures on the warm and cold sides are $T_{+}$ and $T_{-},$ respectively. Typical values may be $T_{+}=30$ ºC and $T_{-}=5$ ºC. It should be noted that the preceding formulas for $p_{\nu }^{+}$ and $p_{\nu }^{-}$ are exact if the flows in all pipe segments are turbulent. There is a small difference if the flow in a segment becomes laminar. This is discussed in the section “Turbulent flow in all pipe segments” (before and after Equation 31(31) $$\mu =0.001\mathrm{(}\text{kg/}\mathrm{(}\mathrm{m}\cdot \text{s))},\hspace{1em}{k}_2=1.571\mathrm{(}\text{kg/}\mathrm{(}\mathrm{m}\cdot \text{s))},\hspace{1em}{D}_{\nu }=0.25\text{(m)}:\hspace{1em}{k}_2\cdot D_{\nu }=0.393$$(31) ).

To elaborate more on the impact of varying viscosities, Figure 11 shows the excess pressure for the network example with five prosumers previously introduced in the section “Bifurcation function for any number of prosumers” and Figure 10. The proportion between the pressures in each node is equivalent to the fourth root factors provided in the third row of Table 1 for $\mu =\mu (5),$ $\mu =\mu (20),$ and $\mu =\mu (35).$

Fig. 11. Excess pressure $p_{\nu }$ at different viscosities for the network example with five prosumers.

Two closed loops, two bifurcation functions

This section investigates the case of network expansion to two closed loops, as shown in Figure 12. There are two connected closed loops involving four parts, each consisting of a sequence of pipes. Parts 1, 2, and 3 make up the first inner closed loop on the right-hand side, while part 2 and the outer part 4 define the second outer loop on the left-hand side. The two loops have part 2 in common. The four parts 1, 2, 3, and 4 are connected to any number $\nu 1,\nu 2,\nu 3,$ and $\nu 4$ of prosumers, respectively.

Fig. 12. Piping grid with two closed loops. The right-hand inner loop consists of three parts, 1, 2, and 3, and the lefthand outer loop of parts 2 and 4.

The sums of prosumer flows $\Sigma 1_{\nu }$ and the mass flows $\dot{m}{1}_{\nu }$ for part 1 become, in accordance with Equation 17(17) $$\boxed{{\dot{m}}_{\nu }={\dot{m}}_0-{\Sigma }_{\nu },\hspace{1em}\nu =0,1,\dots {\nu }_{\mathit{\text{ps}}};\hspace{1em}{\Sigma }_0=0,\hspace{1em}\nu \ge 1:{\Sigma }_{\nu }=\sum _{j=1}^{\nu }{\dot{M}}_j}$$(17) , (39) $$\dot{m}{1}_{\nu }=\dot{m}{1}_0-\Sigma 1_{\nu },\hspace{1em}\nu =0,1,\dots \nu 1;\hspace{1em}\Sigma 1_0=0,\hspace{1em}\nu \ge 1:\Sigma 1_{\nu }=\sum _{j=1}^{\nu }\dot{M}{1}_j$$(39)

There are similar relations for the other three parts: (40) $$\begin{array}{c}\dot{m}{2}_{\nu }=\dot{m}{2}_0-\Sigma 2_{\nu },\hspace{1em}\dot{m}{3}_{\nu }=\dot{m}{3}_0-\Sigma 3_{\nu },\hspace{0.17em}\dot{m}{4}_{\nu }=\dot{m}{4}_0-\Sigma 4_{\nu };\\ \Sigma 2_{\nu }=\sum _{j=1}^{\nu }\dot{M}{2}_j,\hspace{1em}\Sigma 3_{\nu }=\sum _{j=1}^{\nu }\dot{M}{3}_j,\hspace{1em}\Sigma 4_{\nu }=\sum _{j=1}^{\nu }\dot{M}{4}_j,\hspace{1em}\Sigma 2_0=\Sigma 3_0=\Sigma 4_0=0\end{array}$$(40)

The two preceding relations determine the mass flows in all pipe segments with the mass flows in the first segments ($\dot{m}{1}_0,$ $\dot{m}{2}_0,$ $\dot{m}{3}_0$ and $\dot{m}{4}_0$) as unknowns.

There is a mass balance equation at the node where parts 1, 2, and 4 meet: (41) $$\dot{m}{4}_0+\dot{m}{2}_0=\dot{m}{1}_{\nu 1}=\dot{m}{1}_0-\Sigma 1_{\nu 1}\hspace{1em}\Rightarrow \hspace{1em}\dot{m}{2}_0=\dot{m}{1}_0-\dot{m}{4}_0-\Sigma 1_{\nu 1}$$(41)

There is another mass balance equation at the node where parts 3, 2, and 4 meet: (42) $$\dot{m}{3}_0=\dot{m}{2}_{\nu 2}+\dot{m}{4}_{\nu 4}=\dot{m}{2}_0-\Sigma 2_{\nu 2}+\dot{m}{4}_0-\Sigma 4_{\nu 4}=\dot{m}{1}_0-\Sigma 1_{\nu 1}-\Sigma 2_{\nu 2}-\Sigma 4_{\nu 4}$$(42)

The mass balance at the node connected to the accumulator becomes (43) $$\dot{m}{1}_0-\dot{m}{3}_{\nu 3}=\dot{m}{1}_0-\dot{m}{3}_0+\Sigma 3_{\nu 3}=\Sigma 1_{\nu 1}+\Sigma 2_{\nu 2}+\Sigma 3_{\nu 3}+\Sigma 4_{\nu 4}=\Sigma {\dot{M}}_{\mathit{\text{all}}}$$(43)

The mass flow to the accumulator is equal to the sum of all prosumer flows.

Equations 39–42 determine the mass flows in all pipe segments except for two unknowns, $X=\dot{m}{1}_0$ and $Y=\dot{m}{4}_0,$ associated with the two bifurcation loops: (44) $$X=\dot{m}{1}_0,\hspace{1em}Y=\dot{m}{4}_0,\hspace{1em}\dot{m}{2}_0=X-Y-\Sigma 1_{\nu 1},\hspace{1em}\dot{m}{3}_0=X-\Sigma 1_{\nu 1}-\Sigma 2_{\nu 2}-\Sigma 4_{\nu 4}$$(44)

The total pressure change over all segments of the inner loop is a function of the unknown flows $X$ and $Y:$ (45) $$\begin{array}{c}\Delta p_{\text{loop}1+2+3}(X,Y)=\Delta p_{\text{part}1}(X,Y)+\Delta p_{\text{part}2}(X,Y)+\Delta p_{\text{part}3}(X,Y)=\Delta p_{\text{bf\_in}}(X,Y)\\ \end{array}$$(45)

Applying Equation 19(19) $$p_{\nu +1}-p_{\nu }=k_0\cdot \frac{L_{\nu }}{D_{\nu }^3}\cdot f_p\left(\frac{{\dot{m}}_0-{\Sigma }_{\nu }}{k_2{D}_{\nu }}\right),\hspace{1em}\nu =0,1,\dots {\nu }_{\text{ps}};\hspace{1em}{p}_{{\nu }_{\text{ps}}+1}=p_0$$(19) to parts 1, 2, and 3 with mass flows $\dot{m}{1}_{\nu },$ $\dot{m}{2}_{\nu },$ and $\dot{m}{3}_{\nu }$ from Equations 39–44 gives the inner bifurcation function: (46) $$\begin{array}{c}\Delta p_{\text{bf\_in}}(X,Y)=k_0\cdot \sum _{\nu =0}^{\nu 1}\frac{L{1}_{\nu }}{D{1}_{\nu }^3}\cdot f_p\left(\frac{X-\Sigma 1_{\nu }}{k_2\cdot D{1}_{\nu }}\right)+k_0\cdot \sum _{\nu =0}^{\nu 2}\frac{L{2}_{\nu }}{D{2}_{\nu }^3}\cdot f_p\left(\frac{X-Y-\Sigma 1_{\nu 1}-\Sigma 2_{\nu }}{k_2\cdot D{2}_{\nu }}\right)+\\ \\ \hspace{1em}{k}_0\cdot \sum _{\nu =0}^{\nu 3}\frac{L{3}_{\nu }}{D{3}_{\nu }^3}\cdot f_p\left(\frac{X-\Sigma 1_{\nu 1}-\Sigma 2_{\nu 2}-\Sigma 4_{\nu A}-\Sigma 3_{\nu }}{k_2\cdot D{3}_{\nu }}\right)\end{array}$$(46)

The total pressure change over all segments of the outer loop becomes (47) $$\begin{array}{c}\Delta p_{\text{loop}2+4}(X,Y)=\Delta p_{\text{part}2}(X,Y)+\Delta p_{\text{part}4}(X,Y)=\Delta p_{\text{bf\_out}}(X,Y)\\ \end{array}$$(47) (48) $$\Delta p_{\text{bf\_out}}(X,Y)=k_0\cdot \sum _{\nu =0}^{\nu 2}\frac{L{2}_{\nu }}{D{2}_{\nu }^3}\cdot f_p\left(\frac{X-Y-\Sigma 1_{\nu 1}-\Sigma 2_{\nu }}{k_2\cdot D{2}_{\nu }}\right)-k_0\cdot \sum _{\nu =0}^{\nu 4}\frac{L\mathrm{ }\mathrm{ }{4}_{\nu }}{D{4}_{\nu }^3}\cdot f_p\left(\frac{Y-\Sigma 4_{\nu }}{k_2\cdot D{4}_{\nu }}\right)$$(48)

In the second sum. $-\dot{m}{4}_{\nu }=-(Y-\Sigma 4_{\nu })$ shall be used in the argument of $f_p$ in accordance with the sign definitions used in Figure 12. Note also that $f_p(x)$ is an odd function.

The total pressure change over all segments of any closed loop must be zero since the pressure at the start and end refers to the same node. The values of $X$ and $Y,$ where both bifurcation functions are zero, determine the unknown mass flows $\dot{m}{1}_0$ and $\dot{m}{4}_0:$ (49) $$\Delta p_{\text{bf\_in}}(X,Y)=0,\hspace{1em}\Delta p_{\text{bf\_out}}(X,Y)=0\hspace{1em}\iff \hspace{1em}X=\dot{m}{1}_0,\hspace{1em}Y=\dot{m}{4}_0$$(49)

The two unknown mass flows are readily determined by applying a “root” solver two times.

A solver to find the root of a function $F(x)$ is defined by the relations (50) $$F(x)=0\hspace{1em}\iff \hspace{1em}{x}_{\text{root}}=\text{root}(F(x),x)\hspace{1em}\iff \hspace{1em}F(x_{\text{root}})=0$$(50)

The second bifurcation function ($\Delta p_{\text{bf\_out}}$) is zero along a certain “isobar” where $X$ is a function of $Y:$ (51) $$\Delta p_{\text{bf\_out}}(X,Y)=0\hspace{1em}\iff \hspace{1em}{X}_{\text{out}}(Y)=\text{root}\left[\Delta p_{\text{bf\_out}}(X,Y),X\right]$$(51)

This function, $X=X_{\text{out}}(Y),$ is inserted in the first bifurcation function ($\Delta p_{\text{bf\_in}}$). The root (in $Y$) gives $Y_{\text{root}}:$ (52) $$\Delta p_{\text{bf\_in}}(X_{\text{out}}(Y),Y)=0\hspace{1em}\iff \hspace{1em}{Y}_{\text{root}}=\text{root}\left\{\Delta p_{\text{bf\_in}}\left[X_{\text{out}}(Y),Y\right],Y\right\}=\dot{m}{4}_0$$(52)

The other unknown becomes (53) $$X_{\text{root}}=X_{\text{out}}(Y_{\text{root}})=\dot{m}{1}_0$$(53)

There is a close analogy between a pipe network with its flow of water in pipe segments and an electric network with its electric currents. Kirchhoff’s two laws of an electric network state that the sum of the electric currents to a node or junction is zero and that the sum of the electric potentials around any closed loop of electrical resistances is zero.

Discussion and conclusions

The mathematical content of the preceding analyses and solutions is summarized and discussed in this section. In the subsection “Summary of the solution for a single bifurcation loop” the input and all formulas for the solution for a network with a single closed loop of pipes are presented. The flows in the pipes are in our applications normally turbulent. This leads to some simplifications, as discussed in the subsection “Solution for turbulent flows in all pipes with two viscosities.” The case of network expansion to two closed loops shown in the section “Two closed loops, two bifurcation functions” is outlined in the subsection “The case of two bifurcation loops.” The modification to include pipe roughness is presented in the subsection “Impact of pipe roughness.” The model limitations and validations are discussed in the subsection “Model limitations and validation.” Finally, conclusions are outlined in the concluding remarks.

Summary of the solution for a single bifurcation loop

The considered type of network with a warm and cold loop is shown in Figure 1 for three prosumers and in Figure 6 for any number ${\nu }_{\text{ps}}$ of prosumers. The input data required for the full solution are (54) $$\rho ,\mu ,{\mathrm{\text{Re}}}_{\mathit{\text{tr}}}(=2000),{\nu }_{\text{ps}};\hspace{1em}{\dot{M}}_{\nu }\text{for}\nu =1,..{\nu }_{\text{ps}};\hspace{1em}{L}_{\nu },D_{\nu }\text{for}\nu =0,1,..{\nu }_{\text{ps}}$$(54)

Three constants are used, from Equation 13(13) $$k_0=\frac{32{\mu }^2}{\rho }\cdot {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}\text{(N)},\hspace{1em}{k}_1=\frac{0.079}{16}\cdot {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}^{0.75}\mathrm{(}‐\mathrm{)},\hspace{1em}{k}_2=\frac{\pi \mu {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}}{4}\mathrm{(}\text{kg/}\mathrm{(}\mathrm{m}\cdot \text{s))}$$(13) : (55) $$k_0=\frac{32{\mu }^2}{\rho }\cdot {\mathrm{\text{Re}}}_{\mathit{\text{tr}}},\hspace{1em}{k}_1=\frac{0.079}{16}\cdot {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}^{0.75},\hspace{1em}{k}_2=\frac{\pi \mu {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}}{4}$$(55)

The accumulative prosumer flows are also used: (56) $${\Sigma }_0=0,\hspace{1em}\nu =1,\dots {\nu }_{\mathit{\text{ps}}}:{\Sigma }_{\nu }=\sum _{j=1}^{\nu }{\dot{M}}_j$$(56)

The flow-pressure equation for a pipe segment is defined by Equations 3(3) $$\dot{m}=\rho \cdot \pi \cdot {\left(\frac{D}{2}\right)}^2\cdot u\mathrm{(}\text{kg/s}\mathrm{)},\hspace{1em}\mathrm{\text{Re}}=\frac{\rho Du}{\mu },\hspace{1em}{\mathrm{\text{Re}}}_{\mathit{\text{tr}}}=2000$$(3) , 9(9) $$x=\frac{\mathrm{\text{Re}}}{{\mathrm{\text{Re}}}_{\mathit{\text{tr}}}}(‐),\hspace{1em}{k}_0=\frac{32{\mu }^2}{\rho }\cdot {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}\text{(N)},\hspace{1em}{k}_1=\frac{0.079}{16}\cdot {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}^{0.75}\mathrm{(}‐\mathrm{)}$$(9) , 10(10) $$\Delta p=k_0\cdot \frac{L}{D^3}\cdot f_p(x),\hspace{1em}{f}_p(x)=\{\begin{array}{cc}x\hspace{1em}& \left|x\right|\le 1\\ k_1\cdot x\cdot {\left|x\right|}^{0.75}& \left|x\right|>1\end{array},\hspace{1em}x=\frac{\mathrm{\text{Re}}}{{\mathrm{\text{Re}}}_{\mathit{\text{tr}}}}=\frac{\rho D\mathrm{ }u}{\mu {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}}$$(10) , and 12(12) $$\dot{m}=\frac{\pi \rho D^2\cdot u}{4}=\frac{\pi \mu D\cdot \mathrm{\text{Re}}}{4},\hspace{1em}x=\frac{\mathrm{\text{Re}}}{{\mathrm{\text{Re}}}_{\mathit{\text{tr}}}}=\frac{4\dot{m}}{\pi \mu D\cdot {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}}=\frac{\dot{m}}{k_2\cdot D};\hspace{1em}{k}_2=\frac{\pi \mu {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}}{4}$$(12) : (57) $$\Delta p=k_0\cdot \frac{L}{D^3}\cdot f_p(x),\hspace{1em}{f}_p(x)=\{\begin{array}{ll}x\hspace{1em}& \left|x\right|\le 1\\ k_1\cdot x\cdot {\left|x\right|}^{0.75}& \left|x\right|>1\end{array}$$(57)

The bifurcation function gives the pressure drop around the closed loop of pipes as a function of the unknown $\dot{m}={\dot{m}}_0:$ (58) $$\Delta p_{\text{bf}}(\dot{m})=k_0\cdot \sum _{\nu =0}^{{\nu }_{\text{ps}}}\frac{L_{\nu }}{D_{\nu }^3}\cdot f_p\left(\frac{\dot{m}-{\Sigma }_{\nu }}{k_2{D}_{\nu }}\right)$$(58)

The unknown flow in the first pipe segment $\nu =0,$ ${\dot{m}}_0,$ is the value where the bifurcation function becomes zero: (59) $$\Delta p_{\text{bf}}({\dot{m}}_0)=0,\hspace{1em}\left(\sum _{\nu =0}^{{\nu }_{\text{ps}}}\left(p_{\nu +1}-p_{\nu }\right)=p_{{\nu }_{\mathit{\text{ps}}}+1}-p_0=0\right)$$(59)

All pipe mass flows are now determined from the root ${\dot{m}}_0$ in Equation 59(59) $$\Delta p_{\text{bf}}({\dot{m}}_0)=0,\hspace{1em}\left(\sum _{\nu =0}^{{\nu }_{\text{ps}}}\left(p_{\nu +1}-p_{\nu }\right)=p_{{\nu }_{\mathit{\text{ps}}}+1}-p_0=0\right)$$(59) and part of Equation 17(17) $$\boxed{{\dot{m}}_{\nu }={\dot{m}}_0-{\Sigma }_{\nu },\hspace{1em}\nu =0,1,\dots {\nu }_{\mathit{\text{ps}}};\hspace{1em}{\Sigma }_0=0,\hspace{1em}\nu \ge 1:{\Sigma }_{\nu }=\sum _{j=1}^{\nu }{\dot{M}}_j}$$(17) , left: (60) $$\nu =0,1,\dots {\nu }_{\mathit{\text{ps}}}:\hspace{1em}{\dot{m}}_{\nu }={\dot{m}}_0-{\Sigma }_{\nu }$$(60)

The pressures at the nodes are obtained from the recursive relations of Equation 19(19) $$p_{\nu +1}-p_{\nu }=k_0\cdot \frac{L_{\nu }}{D_{\nu }^3}\cdot f_p\left(\frac{{\dot{m}}_0-{\Sigma }_{\nu }}{k_2{D}_{\nu }}\right),\hspace{1em}\nu =0,1,\dots {\nu }_{\text{ps}};\hspace{1em}{p}_{{\nu }_{\text{ps}}+1}=p_0$$(19) : (61) $$p_0=0,\hspace{1em}\nu =0,1,\dots {\nu }_{\text{ps}}:\hspace{1em}{p}_{\nu +1}=p_{\nu }+k_0\cdot \frac{L_{\nu }}{D_{\nu }^3}\cdot f_p\left(\frac{{\dot{m}}_{\nu }}{k_2{D}_{\nu }}\right),\hspace{1em}\left(p_{{\nu }_{\text{ps}}+1}=p_0\right)$$(61)

The excess pressure at the accumulator node, $p_0,$ is zero.

Solution for turbulent flows in all pipes with two viscosities

The equation to determine ${\dot{m}}_0$ for turbulent flow in all pipes is, from Equation 30(30) $$\boxed{\sum _{\nu =0}^{{\nu }_{\text{ps}}}\frac{L_{\nu }}{D_{\nu }^{4.75}}\cdot \left({\dot{m}}_0-{\Sigma }_{\nu }\right)\cdot {\left(|{\dot{m}}_0-{\Sigma }_{\nu }|\right)}^{0.75}=0,{\Sigma }_0=0,\nu \ge 1:{\Sigma }_{\nu }=\sum _{j=1}^{\nu }{\dot{M}}_j}$$(30) , (62) $$\boxed{\sum _{\nu =0}^{{\nu }_{\text{ps}}}\frac{L_{\nu }}{D_{\nu }^{4.75}}\cdot \left({\dot{m}}_0-{\Sigma }_{\nu }\right)\cdot {\left(|{\dot{m}}_0-{\Sigma }_{\nu }|\right)}^{0.75}=0,\hspace{0.17em}{\Sigma }_0=0,\hspace{0.17em}\nu \ge 1:{\Sigma }_{\nu }=\sum _{j=1}^{\nu }{\dot{M}}_j}$$(62)

The mass flow ${\dot{m}}_0$ is independent of the viscosity. The solution for a single prosumer becomes, from Equation 33(33) $${\left(\frac{L_0}{L_1}\cdot \frac{D_1^{4.75}}{D_0^{4.75}}\right)}^{\frac{4}{7}}\cdot {\dot{m}}_0={\dot{M}}_1-{\dot{m}}_0,\hspace{1em}4.75=\frac{19}{4}\hspace{1em}\Rightarrow \hspace{1em}{\dot{m}}_0=\frac{{\dot{M}}_1}{1+{\left(L_0/L_1\right)}^{\frac{4}{7}}\cdot {\left(D_1/D_0\right)}^{\frac{19}{7}}}$$(33) , (63) $${\nu }_{\text{ps}}=1:\hspace{1em}{\dot{m}}_0=\frac{{\dot{M}}_1}{1+{\left(L_0/L_1\right)}^{\frac{4}{7}}\cdot {\left(D_1/D_0\right)}^{\frac{19}{7}}}$$(63)

The pressures at all nodes are given by Equation 35(35) $$\frac{k_0\cdot k_1}{{\left(k_2\right)}^{1.75}}=\frac{32{\mu }^2}{\rho }\cdot {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}\cdot \frac{0.079}{16}\cdot {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}^{0.75}\cdot {\left(\frac{4}{\pi \mu {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}}\right)}^{1.75}=k_{\mu }\cdot \sqrt[4]{\mu },\hspace{1em}{k}_{\mu }=\frac{2\cdot 0.079}{\rho }\cdot {\left(\frac{4}{\pi }\right)}^{1.75}=\frac{0.241}{\rho }$$(35) at the right, and Equations 36(36) $${\tilde{p}}_{\nu }=\frac{p_{\nu }}{\sqrt[4]{\mu }}:\hspace{1em}{\tilde{p}}_{\nu +1}-{\tilde{p}}_{\nu }=k_{.\mu }\cdot \frac{L_{\nu }}{D_{\nu }^{4.75}}\cdot \left({\dot{m}}_0-{\Sigma }_{\nu }\right)\cdot {|{\dot{m}}_0-{\Sigma }_{\nu }|}^{0.75},\hspace{1em}\nu =0,1,\dots {\nu }_{\text{ps}};\hspace{1em}{\tilde{p}}_0=0$$(36) and 37(37) $$p_{\nu }=\sqrt[4]{\mu }\cdot {\tilde{p}}_{\nu },\hspace{1em}{p}_{\nu }^{+}=\sqrt[4]{\mu (T_{+})}\cdot {\tilde{p}}_{\nu },\hspace{1em}{p}_{\nu }^{-}=-\sqrt[4]{\mu (T_{-})}\cdot {\tilde{p}}_{\nu },\hspace{1em}\nu =0,1,\dots {\nu }_{\text{ps}}$$(37) : (64) $$\begin{array}{l}{\tilde{p}}_{\nu +1}-{\tilde{p}}_{\nu }=k_{.\mu }\cdot \frac{L_{\nu }}{D_{\nu }^{4.75}}\cdot \left({\dot{m}}_0-{\Sigma }_{\nu }\right)\cdot {\left|{\dot{m}}_0-{\Sigma }_{\nu }\right|}^{0.75},\hspace{1em}{\tilde{p}}_0=0;\hspace{1em}{k}_{\mu }=\frac{0.241}{\rho };\\ \\ p_{\nu }^{+}=\sqrt[4]{\mu (T_{+})}\cdot {\tilde{p}}_{\nu },\hspace{1em}{p}_{\nu }^{-}=-\sqrt[4]{\mu (T_{-})}\cdot {\tilde{p}}_{\nu },\hspace{1em}\nu =0,1,\dots {\nu }_{\text{ps}}\end{array}$$(64)

It may be noted again that in applications considered in this study, the (nonzero) prosumer flows normally correspond to turbulent flow: $\left|{\dot{M}}_{\nu }\right|>k_2\cdot D_{\nu }$ for all $\nu .$ The preceding equation can then be used to determine ${\dot{m}}_0,p_{\nu }^{+}\text{and}{p}_{\nu }^{-}.$

The case of two bifurcation loops

An example of a pipe network with two closed loops of pipe segments is presented in the section “Two closed loops, two bifurcation functions.” There are two coupled closed loops involving four parts, each consisting of a sequence of pipes. Each one of the four parts is connected to a number of prosumers.

Equation 17(17) $$\boxed{{\dot{m}}_{\nu }={\dot{m}}_0-{\Sigma }_{\nu },\hspace{1em}\nu =0,1,\dots {\nu }_{\mathit{\text{ps}}};\hspace{1em}{\Sigma }_0=0,\hspace{1em}\nu \ge 1:{\Sigma }_{\nu }=\sum _{j=1}^{\nu }{\dot{M}}_j}$$(17) is applicable for all four parts. All mass flows are obtained except the first one (for $\nu =0$). There remain four unknowns to determine. The second closed loop is connected to the first one at two nodes, each with a straightforward mass balance. This leaves two unknown flows $X$ and $Y.$ Each one of the two closed loops is associated with a bifurcation function of $X$ and $Y$ with one sum of pressure differences for each involved part of the loop. The two bifurcation functions are both zero for a certain $X$ and a certain $Y.$ A suitable use of a root solver applied two times determines the two unknowns.

Impact of pipe roughness

Several equations can be adopted to account for the pipe roughness $e$ such as given by Chen (1979), Churchill (1977), Haaland (1983), and Swamee and Jain (1976). In this section, the flow-pressure function in Equation 10(10) $$\Delta p=k_0\cdot \frac{L}{D^3}\cdot f_p(x),\hspace{1em}{f}_p(x)=\{\begin{array}{cc}x\hspace{1em}& \left|x\right|\le 1\\ k_1\cdot x\cdot {\left|x\right|}^{0.75}& \left|x\right|>1\end{array},\hspace{1em}x=\frac{\mathrm{\text{Re}}}{{\mathrm{\text{Re}}}_{\mathit{\text{tr}}}}=\frac{\rho D\mathrm{ }u}{\mu {\mathrm{\text{Re}}}_{\mathit{\text{tr}}}}$$(10) is modified based on the Swamee–Jain equation and the following combined formula is obtained: (65) $$\Delta p=k_0\cdot \frac{L}{D^3}\cdot f_{p_{\text{SW}}}(x,é),\hspace{1em}{f}_{p_{\text{SW}}}(x,é)=\{\begin{array}{cc}x\hspace{1em}& \left|x\right|<1\\ \frac{{\mathrm{\text{Re}}}_{\text{tr}}\cdot x\cdot \left|x\right|}{256{\left(\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}\left(y\right)\right)}^2}& \left|x\right|\ge 1\end{array},\hspace{1em}y=\frac{é}{3.7}+\frac{5.74}{{\left({\mathrm{\text{Re}}}_{\text{tr}}\cdot \left|x\right|\right)}^{0.9}},\hspace{1em}é=\frac{e}{D}$$(65)

The dimensionless flow-pressure function based on the Swamee–Jain equation, $f_{p_{\text{SW}}}(x,é),$ is shown in Figure 13 for different values of relative roughness $é.$ Contrasting Figure 13 with Figure 4 shows that the flow-pressure function for smooth pipes is almost identical. Noticeable differences only occur for very rough pipes, which are not commonly used in applications considered in this study.

Fig. 13. Flow-pressure function $f_{p_{\text{SW}}}(x,é)$ based on the Swamee–Jain equation up to larger x.

Model limitations and validation

The application of the developed method has to be seen in light of some limitations. The first is that the method only considers the main piping loop without dealing much with the prosumers’ components such as service pipes, circulation pumps, heat exchangers, and substation friction losses. The second limitation is related to the representation of the network that must be a single closed piping loop with the same start and end points. Another limitation is that the method currently solves the closed piping loop for a common viscosity along the network length regardless of the different temperature levels at each prosumer. The validation of the newly developed method against three established numerical methods is presented in detail in another article (Abugabbara et al. 2024). In summary, nine validation test cases have been described to simulate different operating scenarios, including, for example, nominal operating conditions, peak summer and peak winter, and network expansion to two closed loops. Results from all modeled cases showed excellent agreement between the new analytical and numerical solution methods for all considered validation metrics.

Concluding remarks

New generations of district heating and cooling systems enable the integration of renewable energy sources and multiple heat prosumers. Due to fluctuations in the heat sources and load variability, each connection point can increase or decrease the network pressure. It is therefore important to ensure safe and adequate network operation by modeling network hydraulics accurately yet quickly.

A new fast analytical method is presented in this article. The method mainly relies on the definition of a simple flow-pressure function used for calculating pressure drops in circular pipes in different flow regimes. Additionally, a bifurcation function is used to determine the unknown mass flow to satisfy network mass balance. The presented solution could be easily modified to include pipe roughness and a solution based on the Swamee–Jain equation for pipe friction factor is provided.

The solution sequence in the newly developed method is straightforward. Only input data with regard to medium properties (density and viscosity), number of prosumers and their corresponding flows, and pipe geometry (length and hydraulic diameter) need to be defined. The bifurcation function and the required constants in the calculations are then defined. Results are finally provided in terms of mass flows in each pipe segment in the loop and the corresponding pressure drops, in addition to the absolute pressure in the nodes.

To the authors’ best knowledge, the presented analytical method for modeling network hydraulics in bidirectional district systems with any number of connected prosumers is the only robust and fast formulation published to date. The developed method can be implemented to account for steady-state or dynamic situations according to the aims of the modeling tool. The method is also capable of modeling meshed network topology consisting of multiple closed loops, which can be realized in future expansion of existing networks. Compared to other established numerical methods, the new analytical method is more favorable due to its explicit nature, ease of implementation, and quick modeling time.

Given the expected increase in the share of pumping energy, a possible future research direction could investigate different optimization techniques to lower the electric energy consumption in DHC systems by implementing the newly developed method.

Nomenclature Acronyms
BLTN	=	bidirectional low temperature network
DC	=	district cooling
DH	=	district heating
DHC	=	district heating and cooling
RES	=	renewable energy sources

Variables
$D$	=	hydraulic diameter [m]
$c_f$	=	friction factor [–]
$e$	=	absolute pipe roughness [m]
$é$	=	relative roughness, $é=e/D$ [–]
$f_p(x)$	=	flow-pressure function based on the Blasius equation [–]
$f_{p_{\text{SW}}}(x,é)$	=	flow-pressure function based on the Swamee–Jain equation [–]
$k_0$	=	constant used for pressure drop calculations [N]
$k_1$	=	constant used for pressure drop calculations [–]
$k_2$	=	constant used for pressure drop calculations [kg/(m·s)]
$L$	=	length [m]
$p$	=	pressure [Pa]
$\dot{m}$	=	mass flow rate in a pipe segment [kg/s]
$\dot{M}$	=	prosumer mass flow rate [kg/s]
$\mathrm{\text{Re}}$	=	Reynolds number [–]
$u$	=	fluid velocity, [m/s]
$x$	=	normalized Reynolds number [–]

Greek letters
$\Delta p$	=	pressure drop [Pa]
$\mu$	=	fluid dynamic viscosity [Pa·s] or [kg/(m·s)]
${\nu }_{\text{ps}}$	=	number of prosumers [–]
$\rho$	=	fluid density [kg/m3]

Subscripts
p	=	pressure
ps	=	prosumer
tr	=	transition to turbulent
$\nu$	=	index for pipe segments and prosumers

Author contributions

Johan Claesson: Conceptualization (equal), Methodology (lead), Software (lead), Visualization (equal), Validation (support), Formal analysis (lead), Investigation (support), Resources (lead), Data curation (equal), Writing – original draft (lead), Writing – review and editing (equal). Jonas Lindhe: Conceptualization (equal), Methodology (support), Software (support), Validation (support), Formal analysis (support), Investigation (support), Resources (support), Data curation (support), Writing – original draft (support), Writing – review and editing (equal). Marwan Abugabbara: Conceptualization (equal), Methodology (support), Software (support), Validation (lead), Visualization (equal), Formal analysis (support), Investigation (lead), Resources (support), Data curation (equal), Writing – original draft (lead), Writing – review and editing (equal). Dennis Johansson: Conceptualization (equal), Resources (lead), Writing – review and editing (equal), Supervision (support), Project administration (equal), Funding acquisition (equal). Saqib Javed: Conceptualization (equal), Investigation (support), Resources (lead), Writing – review and editing (equal), Supervision (lead), Project administration (equal), Funding acquisition (equal).

Disclosure statement

No potential conflict of interest was reported by the authors.