Thermal Transients in District Heating Systems

TL;DR

This paper analyzes the dynamics of heat flux control in district heating systems, focusing on steady velocity flows and their transient behaviors over time, with implications for efficient system operation.

Contribution

It provides a detailed classification of physical phenomena in heat transient propagation and introduces an efficient method to model temperature and heat flux dynamics in networks.

Findings

Transient heat fluxes last tens of minutes to hours.

Turbulent spread of the heat front is a key phenomenon.

Network modeling can be simplified for practical analysis.

Abstract

Heat fluxes in a district heating pipeline systems need to be controlled on the scale from minutes to an hour to adjust to evolving demand. There are two principal ways to control the heat flux - keep temperature fixed but adjust velocity of the carrier (typically water) or keep the velocity flow steady but then adjust temperature at the heat producing source (heat plant). We study the latter scenario, commonly used for operations in Russia and Nordic countries, and analyze dynamics of the heat front as it propagates through the system. Steady velocity flows in the district heating pipelines are typically turbulent and incompressible. Changes in the heat, on either consumption or production sides, lead to slow transients which last from tens of minutes to hours. We classify relevant physical phenomena in a single pipe, e.g. turbulent spread of the turbulent front. We then explain how to describe dynamics of temperature and heat flux evolution over a network efficiently and illustrate the network solution on a simple example involving one producer and one consumer of heat connected by "hot" and "cold" pipes. We conclude the manuscript motivating future research directions.

Equations

$$\forall a\in \mathcal{V}:\underset{\text{mbox}flowbalance}{\underbrace{{\phi }_a=\sum _{b\sim a}{\varphi }_{ab}}}$$

$$\forall \{a,b\}\in \mathcal{E}:\underset{\text{mbox}pressuredroptoflowrelation}{\underbrace{\frac{p_a-p_b}{L_{ab}}=F({\varphi }_{ab})}},$$

$${\partial }_{t}T+V{\partial }_{x}T-D{\partial }_x^{2}T+\gamma T=0,$$

$$D\sim Vd,$$

$$\gamma \sim {\kappa }^{3/4}{V}^{3/4}{d}^{-5/4}{\nu }^{-1/2}$$

$$\rho ({\partial }_{t}u+(u\nabla )u)-\nabla p=\nu {\nabla }^{2}u=\text{mbox}flowdriving,\nabla u=0$$

$${\partial }_t\theta +(u\nabla )\theta =\kappa {\nabla }^2\theta =\text{mbox}heatdriving.$$

$$r_{\nu }={\nu }^{3/4}{d}^{1/4}{V}^{-3/4}.$$

$$r_d={\kappa }^{1/4}{\nu }^{1/2}{d}^{1/4}{V}^{-3/4}.$$

$$q(t;x)=\frac{Q(t;x)}{\rho c_p\pi (d/2{)}^2}\doteq VT-D{\partial }_{x}T,$$

$$\rho c_p\pi (d/2{)}^2\approx 10^3\frac{kg}{m^3}\cdot 4.2\cdot 10^3\frac{W\cdot s}{kg\cdot K}\cdot \pi (0.25\cdot m{)}^2\approx 8.2\cdot 10^6\frac{W\cdot s}{K\cdot m},$$

$$\forall x\in [0,L]:$$

$$T(0,x)=T_{st}(x)\doteq \frac{2q_0}{V+\sqrt{V^2+4D\gamma }}\exp \left(\frac{V-\sqrt{V^2+4D\gamma }}{2D}x\right).$$

$$\tilde{T}(t,x)=T(t,x)-T_{st}(x).$$

$${\partial }_t\tilde{T}+V{\partial }_x\tilde{T}-D{\partial }_x^2\tilde{T}+\gamma \tilde{T}=q(t)\delta (x).$$

$$\forall t>0:V\tilde{T}(t,0)-D{\partial }_x\tilde{T}(t,0)=q(t),$$

$$\tilde{T}(t,+\infty )=0,$$

$$\forall x\ge 0:\tilde{T}(0,x)=0.$$

$$\tilde{T}(t,x)=\underset{0}{\overset{t}{\int }}d{t}^{\prime }q(t^{\prime })\frac{\exp \left(-\gamma (t-t^{\prime })-\frac{(x-V(t-t^{\prime }){)}^2}{4D(t-t^{\prime })}\right)}{\sqrt{4\pi D(t-t^{\prime })}}.$$

$$\forall \alpha \in \mathcal{E}:V_{\alpha }=\frac{{\varphi }_{\alpha }}{\pi (d_{\alpha }/2{)}^2\rho }$$

$$\forall a\in \mathcal{V}:\forall \alpha \in \text{mbox}out(a):q_{a\to \alpha }=G_{a\to \alpha }(q_{a;in};Q_a).$$

$$\underline{\text{mbox}Branchingnodewithperfectmixing}:$$

$$\underline{\text{mbox}Standardproduction/consumptionnode(singleinput,singleoutput)}:$$

$$\underline{\forall \alpha =(a,b)\in \mathcal{E}}:$$

$$q_{\alpha \to b}=q_{a\to \alpha }\exp \left(\frac{V_{\alpha }-\sqrt{V_{\alpha }^2+4D_{\alpha }{\gamma }_{\alpha }}}{2D_{\alpha }}{L}_{\alpha }\right),$$

$$\underline{\forall \alpha =(a,b)\in \mathcal{E}}:q_{\alpha \to b}(t)=$$

$$\underset{0}{\overset{t}{\int }}d{t}^{\prime }{q}_{a\to \alpha }(t^{\prime })\left(V_{\alpha }+\frac{L_{\alpha }}{2(t-t^{\prime })}\right)\frac{\exp \left(-{\gamma }_{\alpha }(t-t^{\prime })-\frac{(L_{\alpha }-V_{\alpha }(t-t^{\prime }){)}^2}{4D_{\alpha }(t-t^{\prime })}\right)}{\sqrt{4\pi D_{\alpha }(t-t^{\prime })}}.$$

$$T(t,0)=T_{st}(0)+\tilde{T}(t,0)=T^{(in)}(t),$$

$$\hat{\mathcal{L}}(s;x)T_s(x)=T_{st}(x),$$

$$\hat{\mathcal{L}}(s;x)\doteq s+\alpha +V\frac{d}{dx}-D\frac{d^2}{d{x}^2},$$

$$T_s(0)=T_s^{(in)},T_s(\infty )=0,$$

$$T_s(x)\doteq {\int }_0^{\infty }dt{e}^{-ts}T(t,x),$$

$$T_s^{(in)}\doteq {\int }_0^{\infty }dt\exp (-ts)T^{(in)}(t).$$

$$\displaystyle\Delta^{(lin)}(t)=\Delta_{*}\left\{\begin{array}[]{cc}\frac{t}{\tau},&0\geq t\leq\tau\\ 1,&\tau\leq t\end{array}\right.;\quad\Delta^{(lin)}_{s}=\Delta_{*}\frac{1-e^{-\tau s}}{s^{2}\tau};$$

$${\Delta }^{(exp)}(t)={\Delta }_{\ast }\left(1-e^{-\frac{t}{\tau }}\right);{\Delta }_s^{(exp)}=\frac{{\Delta }_{\ast }}{s(s\tau +1)}.$$

$$T_s(x)={\tilde{T}}_s(x)+\frac{T_{st}(x)}{s},$$

Full text

††journal: Energymytitlenotemytitlenotefootnotetext: Thermal Transients in District Heating Systems: Physics Modeling for better Control

url]https://sites.google.com/site/mchertkov/

Highlights

Advection-Diffusion-Loss of heat in district heating networks is analyzed.

• 2.

Parameters of the basic model follow from turbulence phenomenology.

• 3.

Superposition of running and spreading fronts forms a typical transient.

• 4.

We suggest an efficient computational scheme to describe the transients.

1 Introduction

1.1 Smart Systems

Conceptual development of ”smart grids” has started in Power Systems (PS) [1, 2], however extensions of the concept to other energy infrastructures, such as Natural Gas Systems (NGS) [3, 4, 5] and District Heating (and/or cooling) Systems (DHS) [6, 7], are now considered as well. In particular, smart DHS are of a special interest to countries with significant seasonal variations, such as Northern European countries, Russia and potentially USA and Canada. Even countries with milder climate, such as Southern European countries, which do not rely on DHS with a significant (city-scale) footprint, start to reconsider these practices. The city of Torino (Italy) is one impressive example of this type where a very modern system was built during the last 20 years from a scratch [8]. One may expect, based on this and other recent European examples, that adaptation of this emerging technology to other countries, e.g. USA which so far stay immune to these developments, should be expected soon.

Smartness of a DHS should show itself in how the system is monitored (data), described (physics), controlled (through optimization built on data/measurements and aware of physics) and planned (all of the above). Smart extensions, design and/or re-design of DHS should include long term planning (15-30 years), short term planning (months to years) and operational scheduling/planning (hours to days). Variety of planning and operational problems/formulations should account for uncertainties, take advantage of new technologies and devices which are to be installed in the systems, and build a modern operational version of the system capable to sustain disruptions and emergencies. In combination, these smart developments should be beneficial because of savings, efficiency, convenience and energy security they are expected to bring in. Specific to DHS [9, 10] one aims to (a) reduce losses by lowering temperature in the system and slowing flows, (b) take advantage of active consumers (demand response), (c) add new smarter devices (pumps, storage, co-generation), (d) rely extensively on sensing and measurements of the mass and heat flows through out the system, (e) co-optimize thermal, electric and gas infrastructures.

1.2 New challenges

To advance towards completing the tasks outlined (very schematically) above will face many challenges, in particular:

• (1)

Mass transport in the case of an incompressible fluid chosen as the heat carrier (typically water) establishes fast, through a transient propagating with the speed-of-sound, however thermal profile through out the system settles much slower. Main foundations of the steady mass and heat transport in the DHS were established in 80ies and 90ies, see e.g. [11, 12]. More recently the researchers have focused their attention on analysis of thermal transients [13, 14, 15, 16, 17, 18], e.g. addressing the question of how the heat fronts propagate through the system. In these studies the authors investigated how the temperature field, and thus the heat fronts, are advected by the flow (velocity) along the pipes. Different effects were considered. In particular, [13] reported significant dependence of the thermal front spread in response to a rapid change in the flow velocity. It was emphasized in [14] that effects of the systems geometry (especially bends and fittings) influences the thermal front speed and spread rather significantly. Effect of the rapid change of the inlet temperature of the source on the thermal front was analyzed in [15]. Diffusive spread of the front due to turbulence was mentioned and discussed as potentially important in [13, 14, 15], however it was in the end ignored in the simulations as producing a small effect in the relatively small systems (with spatial extent of 1-2km or less). The turbulent spread of the front was also not accounted in the simulations of [16, 17, 18]. In fact this is not surprising as the systems analyzed in [16, 17, 18] were too small for the turbulent diffusion effects to show a significant prominence.

Even thought it is perfectly reasonable to ignore the turbulent diffusion effects in small systems, or systems where the flows are slow/laminar 111 Note in passing that all the aforementioned studies, reported in [13, 14, 15, 16, 17, 18], were actually conducted in a turbulent regime with sufficiently large Re-number., the turbulent effects should certainly be accounted for in larger systems operated in the turbulent regime (of high Re-number). In other words when the system is turbulent (which is the common case in majority of the district heating systems) and sufficiently large the spread of the front becomes significant enough to cause troubles with timely delivery of services (heat) and with providing sufficiently accurate model for optimal dispatch and control of the heat.

• (2)

Models of DHS should be probabilistic, and as such should account for uncertainty and stochasticity in model parameters, modeling of consumers and producers of heat, and also in modeling of coupling to other infrastructures and (even more generally) other degrees of freedom. Some preliminary work on statistical modeling of DHS was reported by one of the authors of this manuscript [19, 16, 20], however many of relevant statistical questions are yet to be posed and discussed. Note that adapting to DHS statistical and uncertainty-handling methods, recently developed for related problems in power systems, see e.g. [21, 22], is of a great interest.

• (3)

Measurements-based monitoring of the system state and model validation are important tasks which were posed and discussed in the literature preliminarily [19]. Further work aimed at developing scalable practical algorithms of both on-line and off-line type resulting in high-quality and reliable predictions is needed.

• (4)

The last, but not the least, construction of a comprehensive control, optimization and planning paradigm, built on methods and techniques brought in from other science and engineering disciplines to resolve problems and issues mentioned above, constitutes an over-reaching hyper-challenge.

In this paper we mainly focus on addressing the primary challenge #1, i.e. developing dynamic models for the heat transport over large systems (of a modern city scale) operating at high Re numbers which are to be used in the future in tests and validations related to challenges #2-4.

1.3 Operational Modeling

Most important characteristics of the heat flow spatio-temporal evolution within DHS is that it is never steady. This is in contrast with the mass flow, which is almost always steady or quasi-steady, under exception of very short transients when it changes and then settles with the speed of sound. It is appropriate to note here that design of the basic principles of the heat flux control has a significant effect on how one models the thermal/heat transients. We distinguish two principally different paradigms. Control-by-temperature which is adopted in Nordic countries and Russia where the first large (city scale) DHS systems were designed a century ago. In these systems mass flow is adjusted rarely, and even then mass flow transients are settled in the matter of seconds and maximum minutes. Main controls are implemented via change of temperature at the heat sources, that is at the major Combined Heat & Power (CHP) plants and big boilers. The temperature set points are changed through out the day to adjust to the change in the heat demand. Change of the temperature set point generates a heat front propagating with speed of the mass flow which varies from tenth of meters per second to a few meters per second. Propagation of the heat front through a system (where it does not diffuse yet and is seen as a clear front) may have a significant spatial extend – from kilometers to tens of kilometers, lasting from minutes to hours. This raises the questions we aim to address in the paper – of controlling delays in heat delivery and quantifying details of the slow transients. As recognized relatively recently modeling delays, and thus thermal transients while accounting not only for ballistic propagation of the heat fronts but also for diffusive spread of the fronts, is significant [13, 14, 15, 16]. Notice that modern systems build in S. European countries have adopted another operational principe – control-by-mass flows when the amount of heat flux is controlled through adjustment of the flow velocity [8]. In such systems flow velocity and heat fluxes are settled (speed of sound) fast. In the heat delivery (hot) part of such systems temperatures are kept constant through out the day, while in the return (cold) part of the system temperature varies depending on the amount of heat consumed by customers. Therefore propagation and diffusive spread of the heat fronts in such mass-flow-controlled system may be of a significance only in the return (cold) part of the system. Note also that it is natural to expect that the two main operational principles will be merged and combined in smart future DHSs.

1.4 Layout of Material and Our Main Results

Layout of material in the manuscript and our main results reported in this manuscript are as follows:

We briefly review standard description of the steady hydro (mass) flow in the pipes averaged over turbulent fluctuations, in Section 2.

• 2.

Basic equation describing dynamics of temperature in the steady mass flow is discussed in Section 3. We first introduce the basic thermo-advection-losses equation and then, assuming that the mass flow of the incompressible carrier (water) in the DHS is turbulent, we describe in Section 3.1 how standard turbulence theory phenomenology can be used to argue validity of the basic equation. We present parametric estimates for the turbulent diffusion and heat losses coefficients, both averaged over spatial (cross-section of the pipe) and temporal (cross-section of the pipe times inverse flow velocity) as a function of mean flow velocity. This analysis results in a conclusion that the 1+1 (space+time)- dimensional advection-diffusion-losses equation, supplemented by initial and boundary conditions, are appropriate for modeling thermal transient in an individual pipe and system of pipes. Quantitative estimates for the parameters correspondent to practical DHS system are given in Section 3.2. We also discuss estimation for (and importance of) the heat flux in DHS in Section 3.3.

• 3.

In Section 4 we solve the advection-diffusion-losses equation, introduced and discussed in Section 3, in the setting of a single pipe. Relevant additional/auxiliary material is also placed in A.

• 4.

We discuss extension of the modeling described in Section 4 for a single pipe to a DHS network in Section 5, first addressing the static case in Section 5.1 and then dynamics/transient case in Section 5.2. The network solution is illustrated in Section 5.3 on example of the network consisted of one producer of heat and one consumer of heat connected by two pipes, hot and cold.

• 5.

Section 6 is reserved for discussions of the results, conclusions and path forward (future tasks).

2 Hydro (mass-) flow equations

Description of the hydro-thermal flow equations is split in two parts. First, we describe equations governing mass flows over the system pipes. Then, we move to description of the temperature field spatio-temporal spread over the system affected by the flow. We will not consider here the speed-of-sound transients which are settled in a matter of seconds (or faster). We will also assume that the heat carrier is a single-phase fluid (water) with the pipes fully filled (no gas component present) - thus ignoring extreme regimes when the water can vaporize or be in the water-vapor (two-phase) regime. Overall, inncompressibility of the single-phase flow with the focus on physical processes, which are slower than the speed-of-sound effects, means that we are in the balanced regime, i.e. in response to consumption change there should be an instantaneous (within the approximation) adjustment of the production. (The situation is quite different in compressible flows, e.g. describing transport of natural gas in the pipes, where the so-called line pack effect, resulting in variability of the total amount of gas contained within the pipes, takes place.) Under these assumptions, the set of balanced quasi-static equations, explaining distribution of the mass-flow and pressure over the entire system (network) of pipes, becomes

[TABLE]

where ${\cal G}=({\cal V},{\cal E})$ is the undirected graph of the pipe network, and $a\in{\cal V}$ and $\{a,b\}\in{\cal E}$ state that node $a$ and edge $\{a,b\}$ are taken from the set of nodes, ${\cal V}$, and set of edges, ${\cal E}$, respectively. Other characteristics entering Eqs. (1,2) are as follows. $\varphi=(\varphi_{a}|a\in{\cal V})$ is the vector of the mass-flow injections/consumptions which is globally balanced (in view of the comments above), $\sum_{a\in{\cal V}}\varphi_{a}=0$. $\phi=(\phi_{ab}=-\phi_{ba}|\{a,b\}\in{\cal E})$ is the vector of mass flows over pipes, where each component is odd with respect of the (edge/pipe) index permutations. Eq. (2) describes relation between the pressure drop at the two ends of the pipe, $p_{a}-p_{b}$, per length of the pipe, $L_{ab}=L_{ba}$, and the steady mass flow along the pipe, $\phi_{ab}$. $F(\phi)$ is a nonlinear function representing effect of turbulent friction. According to the standard fluid mechanics estimation, $F(\phi)\sim\phi|\phi|$, however some other modifications of the simple quadratic dependence [23] may be more appropriate to model (phenomenologically) effects of pumps, flow regulators and other control devices.

3 Temperature (thermal) advection-diffusion equations: the case of a single pipe

Consider a single pipe. Mass flow is related to density, $\rho$, and (averaged over cross-section) velocity, $V$, according to $\phi=\rho V\pi(d/2)^{2}$, where $d$ is the diameter of the pipe (assumed constant along the pipe). Assuming also that thermodynamic variations (of density and temperature) are small, the constancy of $\phi$ along the pipe translates into constancy of the flow velocity along the pipe.

Coarse-grained, and already averaged over cross-section of a pipe dynamics of a thermal field in a steady (but turbulent) flow is controlled by the following advection-diffusion equation

[TABLE]

where $T(t;x)$ is the temperature field, possibly dependent on time, $t$, and position along the pipe, $x$, and counted from an ambient temperature (here “ambient” should be counted not as an outside temperature but temperature of the pipe, under insulation coat); $D$ is the eddy (turbulent) diffusion coefficient; $\gamma$ is the coefficient of the thermal heat loss through the walls of the pipe (assumed kept at the ambient temperature). $D$ and $\gamma$ are related to $V$, as well as to the geometric parameters of the pipe, and in the turbulent regime can be estimated (through straightforward, a-la Richarson-Kolmogorov-Obukhov-Kraichnan (RKOK), phenomenology) as

[TABLE]

where $\kappa$ and $\nu$ are thermal diffusion coefficient and kinematic viscosity of the water, respectively. We provide brief description of the RKOK picture of flow and heat profile leading to Eqs. (1,2) and Eqs. (3) and resulting in the estimations of eddy diffusivity (4) and thermo-losses (5) respectively, in the next Subsection 3.1.

Realistic estimates for the eddy-diffusivity (4) and thermo-loss (5) coefficients in DHS pipes are given in Subsection 3.2.

3.1 Richardson-Kolmogorov-Obukhov-Kraichnan picture of flow and heat profile in a turbulent pipe

Basic microscopic equations governing dynamics of incompressible velocity, $u$, and temperature, $\theta$, at a moment of time $t$ and at a spatial position, $r$, within the pipe (thus in three dimensions) are (schematically)

[TABLE]

The equations should be complemented by respective boundary conditions. Averaging these equations over turbulent fluctuations in the pipe results in the steady mass flow Eqs. (1,2) and coarse-grained thermo-transport Eqs. (3) for the average (over cross-section) velocity, $V$, and temperature, $T$, thus dependent (potentially) only on the position along the pipe.

To estimate eddy-diffusivity and heat loss coefficients entering the coarse-grained description of Eq. (3) one needs to discuss underlying picture of turbulent boundary layers, schematically illustrated/summarized in Fig. (1). There are actually two boundary layers, one related to velocity – so-called viscous boundary layer – the flow is laminar within the layer, and the other called diffusive boundary layer – within which thermal diffusion dominates turbulent diffusion/advection. Widths of the two boundary layers are estimated within the standard Kolmogorov-Obukhov (KO) phenomenology [24] as follows.

Velocity self advection with viscosity and assuming that velocity fluctuations at the scale, $r$, from within the range of scales bounded from above by $d$ and from below by $r_{\nu}$ scales as $\delta v_{r}\sim V(r/d)^{1/3}$) and thus

[TABLE]

We assume that $\nu>\kappa$ (in water the so-called Prandtl number, defined as $Pr=\nu/\kappa$, ranges from $13$ at $0C$ to $\sim 2$ at $100C$.) Then the thermal diffusion boundary layer is found within the larger viscous boundary layer. Diffusion dominates close to the walls, i.e. within the diffusion layer, while advection takes over on the interface with the bulk. Estimation for the width of this layer, $r_{d}$, from the Obukhov-Batchelor-Kraichnan (OBK) theory of scalar transport in turbulence [24] (assuming that velocity fluctuations in the range of scales bounded from above by $r_{\nu}$ and by $r_{d}$ is $\delta v_{r}\sim V(r_{\nu}/R)^{1/3}(r/r_{\nu})$ and balancing turbulent advection and thermo-diffusion term at $r_{d}$) is

[TABLE]

The general picture of the microscopic temperature profile, $\theta$, across the cross-section of the pipe is roughly as follows. $\theta$ is constant in the central portion of the pipe and then it drops to zero in the thermo-diffusive layer, where thermal flux (losses to the wall) balances thermal diffusion. Therefore, balancing integral of the two terms over circumference of the pipe one estimates, that $\gamma d\sim\kappa/r_{d}$, resulting in Eq. (5).

This concludes our brief recount of the KO-OBK phenomenology as applied to turbulent incompressible flows and scalar transport in a pipe. (See [24] and references there in for more details.) In the following Subsection we present estimation for the thermo-flow characteristics typical for pipe flows in DHS.

3.2 Effective Parameters for an Exemplary Pipe

[TABLE]

Parameters corresponding to exemplary cases from Russia [16], and Denmark [13] are shown in the Table 3.2. Here, $L$ is the length of the pipe/system, $d$ is diameter of the pipe, $V$ is the flow velocity, $D\sim Vd$ is the eddy/turbulent diffusion coefficient, $Re=D/\nu$ is the Reynolds number where $\nu\approx 5\cdot 10^{-7}m^{2}/s$ is the kinematic viscosity of the water (at 60C), $\bar{t}=L/V$ is the expected time of propagation (e.g. of the heat front) through the system, $\Delta L=\sqrt{D\bar{t}}=\sqrt{dL}$ uncertainty in the front position due to eddy diffusivity, $\Delta t=\Delta L/V=\sqrt{dL}/V$ uncertainty in the time of propagation through the system due to turbulence, $\gamma$ is turbulent thermal loss coefficient, estimated according to Eq. (5, and one takes $\kappa\approx 1.8\cdot 10^{-7}\ m^{2}/s$ for water (at 60C).

3.3 Heat flux

An important characteristic of interest is the heat flux transferred along the pipe at the position, $x$, and the moment of time, $t$:

[TABLE]

where the two terms on the right hand side of Eq. (10) correspond to the advective contribution (which may also be called mass flow or hydraulic contribution) and (turbulent, eddy-) diffusive or convective contribution; $\rho c_{p}\pi(R/2)^{2}$ sets the units for the heat flux, so that $Q$, representing power, is measured in [Watts], while $q$ is measured in [m K/s]; $\rho$ is the carrier/fluid density, $c_{p}$ is the heat capacity of the carrier/fluid, and $\pi R^{2}/4$ is the area of the pipe cross-section; Eq. (11) gives an exemplary estimate for water. Notice, that the split on advective/hydraulic vs diffusive/conductive terms is formal in the turbulent regime where, according to Eq. (4), the two mechanisms originate from the same physics and cannot be separated. Notice that in the normal operational regime heat flow goes along the mass flow and also against the temperature gradient, i.e. the two terms on the rhs of Eq. (10) work in unison, even though one expects that the first advective/hydraulic term is significantly larger, since $L\gg R$. (In principle the small advective correction may be of the opposite sign in the system where water in the pipes is kept at the temperature lower than the ambient/outside temperature. However, these cooling regimes are exotic, if realistic at all.)

4 Cauchi problem: formulation and solution

In this Section we discuss solution of the thermal transport Eq. (3) describing propagation of a thermal front through a district heating pipe.

Without loss of generality we assume that initially, i.e. at $t=0$, temperature profile of the temperature in a single pipe of length $L$ is the steady time-independent solution of Eq. (3) correspondent to injection into the pipe of the flux $q_{0}$ at the inlet, $x=0$,

[TABLE]

Given similarity to the so-called Shukhov equation 222 Shukhov is Russian engineer-polymath, scientist and architect [25] who has suggested to model temperature decay in steady pipe flows due to heat losses via an exponential function of the distance (along the pipe). General discussion of the Shukhov formula for dependence of the temperature drop on other characteristics and properties of the flow can be found in [26, 27]. See also [28, 29, 30, 31] for extensive discussions, derivations and validations of the Shukhov formula in the regimes of interest for water-based district heating systems. we call Eq. (12) the generalized Shukhov equation. (To derive the original Shukhov equation from the generalized Shukhov equation, one should simply replace $D$ in Eq. (12) by zero.)

Denoting

[TABLE]

deviations from the steady profile driven by the action of the additional (to $q_{0}$) heat source, $q(t)$, switched on at $t=0$ and growing to a constant, one arrives at the following equations

[TABLE]

We recap that this setting corresponds to enforcing for $\tilde{T}$ the flux Boundary Condition (BC) at the heat injection inlet of the pipe 333 Physical consideration described above justifies this choice of the constant flux condition at the origin. Notice that imposing other BC at the origin, e.g. so-called Neumann BC correspondent to maintaining the temperature constant at the inlet (and not the flux) may also be realistic for some, less common, heat injection sources. Even more generally a mixed BC, maintaining a linear combination of temperature and of the heat flux is also a possibility. Comprehensive discussion of how to model the heat source condition goes beyond the scope of the paper. The authors believe that this important question should be resolved in the future through an accurate model calibration against experimental/field data.,

[TABLE]

zero asymptotic condition at $x={+}\infty$,

[TABLE]

and zero initial condition

[TABLE]

Explicit solution of Eq. (14,15,16,17) becomes

[TABLE]

Advancement of the heat front, its diffusive spread and the follow up transient past the front settling into a steady profile are illustrated in Fig. (2). Discussion of other possible choices of the BC (see the footnote), e.g. corespondent to maintaining the temperature constant at the inlet, is discussed in A.

5 Temperature (thermal) advection-diffusion equations in a multi-pipe system (network)

Consider stationary flows and, moreover, ignore dependence of the flows on density and temperature, assuming that changes in the latter are sufficiently small to have a significant thermodynamic effect on density, $\rho$, considered known and constant through out the multi-pipe system at the time scale of consideration (through out few hours). Then, the system of Eqs. (1,2) are separated (from dynamics of the temperature field and heat fluxes) and shall be re-solved first. Output of the hydro-static analysis is the set of velocity values and mass flows associated with each edge (pipe) of the system

[TABLE]

thus found and assumed known. Here in Eq. (19) $d_{\alpha}$ is the diameter of the pipe $\alpha$ and $\rho$ is the density of fluid/water. One also assumes that $\alpha$ marks a label of a directed edge oriented along the direction of the flow of the given pipe and that the diameter of the pipe is constant (does not change along the pipe), and thus the velocity of the incompressible flow is constant (along the pipe) too.

The hydro-static description of flows by Eq. (1,2,19) is to be complemented by the advection-diffusion equations discussed in the following Subsection.

5.1 Formulating and Solving Stationary Advection-Diffusion-Losses Equations over the Network

Consider a thermo-static case when all the temperatures and heat flows are stationary. Temperature field distribution over a pipe for given boundary conditions on temperature (or heat flux) was described above, thus what remains to be done is to complement this “along the pipe” description with providing details on the temperature change at other “nodal” elements of the system.

Layout of an actual/technological system is rather involved, see e.g. Fig. (3) showing the base elements and an exemplary computational scheme. However, aiming to describe only principal modeling we will not attempt to reproduce all the practical details here, instead we will focus on a somehow simplified view of a district heating network containing only heat branching/mixing nodes and producing/consuming nodes, thus skipping to describe such practically important elements of the district heating system as pumps, heating chambers, etc. In the following we will work primarily with the single-pipe representation where hot and cold pipes are considered separately. See e.g. top sub-figures/portions in Fig. 4.

In fact, aiming to have a simple modular description we will not differentiate between branching nodes and active nodes, combining these in a generalized node, $a\in{\cal V}$, each represented by a set of incoming pipes, $\alpha\in\mbox{in}(a)\subseteq{\cal E}$, and outgoing pipes, $\alpha\in\mbox{out}(a)\subseteq{\cal E}$, where ${\cal E}$ represents the set of directed (oriented along the flow) edges/pipes, and ${\cal V}$ represents the set of nodes. Then heat fluxes of the pipes outgoing a generalized node, $a$, are expressed via incoming heat fluxes, $q_{a;in}=(q_{\alpha\to a}|\alpha\in\mbox{in}(a))$, and of the production/consumption of heat, $Q_{a}$, at the node

[TABLE]

The $G$-functions, representing models of the nodes, can be rather different for different nodes. Here we present some most popular exemplary models. A perfect mixing model of the branching node without injection/consumption of heat, i.e. with $Q_{a}=0$, is

[TABLE]

A standard model of a production/consumption node, $a$, with one input pipe, $\mbox{in}(a)=\alpha$ and one output pipe, $\mbox{out}(a)=\beta$ is

[TABLE]

where $Q_{a}$ is heat injected (then $Q_{a}>0$) or consumed (then $Q_{a}<0$) at the node $a$.

According to Eq. (12), which is the stationary solution of the basic Eq. (3), heat flux at the outlet of the pipe, as expressed via the heat flux at the inlet, is

[TABLE]

where (as before) the in- and out- sides of the pipe are defined with respect to direction of the mass flow.

Eqs. (20), where the respective G-functions can be chosen, for example, according to Eqs. (21) and Eqs. (22), supplemented by Eqs. (23) and defined over a single-connected district heating network, set up an example of the steady “network heat flow” problem:

Input:****

• (a)

Heat flows at all the pipes leaving all the production nodes of the network;

• (b)

Powers consumed at all the consumer nodes;

• 2.

Output:****

• (a)

Heat flows at all other locations (not given within description of the input);

• (b)

Temperatures at all network locations;

• (c)

Powers injected at the production nodes.

Consistently with what was discussed in the literature, see e.g. [23], this set up has a unique solution. Indeed, to verify this we need to start from generators and solve the set of Eqs. (20,21,22,23) recursively, one by one, simply following the directed (in the single-pipe representation) graph of the district heating network. Notice that since the system of equations is linear, re-scaling the solution (one which was just found propagating relations along the directed graph) by multiplying all heat flows, power injections/consumptions and temperature by the same positive constant number will generate another valid solution of the system of Eqs. (20,21,22,23). This one-parametric degree of freedom can be useful for finding solutions of the system of equations in the settings where boundary conditions (input characteristics in the description above) arrange differently. For example, it may be more practical to consider a setting where heat flows at all the pipes leaving all the generators of the network are replaced by powers injected at all the generators of the network. (These two different sets of input parameters are in the one-to-one relation.) Another option for a valid input set is to replace some (or all) of the injected/consumed powers at the producers/consumers by temperatures - one per node (say outgoing for producer and incoming for consumer).

5.2 Formulating and Solving Dynamic/Transient Advection-Diffusion-Losses Equations over Network

Generalization of the static version of the ”heat flow” problem posed in the preceding Subsection to the dynamic/transient case is straightforward - one just needs to substitute Eqs. (23) by the following solution of the dynamic equations derived from Eq. (18)

[TABLE]

Let us emphasize that explicit form of Eq. (24) in quadratures makes the dynamic problem as efficient computationally as the static problem.

5.3 Dynamic Network Illustration

We demonstrate the quality and efficiency of our dynamic/transient advection-diffusion-losses modeling on the simple example of one producer and one consumer connected by two pipes (hot and cold) shown in Fig. 4 (right). We choose parameters for the two pipes the same and correspondent to the single pipe illustration shown in Fig. 2. The results for two different heat injection/consumption settings are shown in Figs.  (5,6). In the first case injection and consumption of heat are started simultaneously, while in the second case consumption of heat starts with a delay. Comparison of the two examples illustrates benefits and importance of the dynamic modeling accounting for time delays and details of the heat front spreading as it advances through the network.

6 Conclusions and Path Forward

In this manuscript we started to explore new phenomena related to the temporal dynamics in the District Heating Systems. Specifically, we have made the following set of statements and observations:

Dynamics of the thermal field in DHS pipes, analyzed over distances which are longer than diameter of the pipe and at temporal scales which are longer than the sound-wave transients, is governed by the set of coupled linear advection-diffusion-heat loss equations. Parametric dependence of the turbulent diffusion and loss coefficients on mean velocity of the turbulent flow in a pipe and microscopic characteristics of the heat carrier (water) were estimated, according to standard turbulent phenomenology, in the practical regime of large Re-numbers.

• 2.

Solutions of the advection-diffusion-heat loss equations were analyzed under rather general initial and boundary conditions. Description of a heat front, initiated by a change of temperature or heat flux at the heat source or sink, was in the focus of our analysis which captures principal qualitative features of the phenomena, such as diffusive spread of the front. General solution was presented in quadratures and compared with direct simulations. The quadrature representation of the general solution in the form of an integral convolution of the heat kernel with the function representing the heat source(s) is explicit.

• 3.

We have explained, in principle and on a simple example, how the single pipe results allow efficient and accurate computations of thermal transients in a network. The problem reduces to explicit recurrency/propagation of solutions from the heat source to consumers and back thus bypassing solving the system of computationally expensive Partial Differential Equations (PDEs) over the network.

In relatively short terms the results presented in the manuscript call for further research and exploration of the following topics:

Demonstration of the power of the approach on a realistic district heating network. This should also include calibration of the freedom in modeling heat sources/producers and sinks/consumers.

• 2.

Extension of the network-wide computations to variety of the operational regimes of interest. In this manuscript we have focused solely on describing regime when the mass/velocity flows are kept steady. However, the methodology is extendable to other regimes of interest, e.g. where the heat flux is regulated through control/change of the mass flow.

• 3.

Accounting for statistical effects and uncertainty.

Further down the road we plan to utilize the approach and address more complex problems, such as

Data driven validation of the network parameters and robust modeling. District heating systems are subject to relatively fast changes acquired through operational wear and tear and related maintenance services. This is typically seen through a rather significant change of the key elements of pipes’ and devices’ parameters acquired in the time frame of a month (especially during the active/cold season). Thus, effective diffusion coefficients and effective heat-loss coefficients, influencing the thermal dynamics in the network, will change on this time scale and may be uncertain. There are two approaches we plan to develop to account for the complications. First of all, we will be developing efficient data driven approaches which allow to reconstruct parameters from the minimal amount of field measurements. Second, even if the data driven reconstruction is employed one still do not expect to reduce the parameter uncertainty completely. Then we will rely on the probabilistic techniques to describe how uncertainties within the model parameters affects description of the heat transfer through the network.

• 2.

Development of controls for improved efficiency and security. Scalability of the thermal transient computations which we started to develop in this manuscript should allow to explore variety of available controls efficiently. In addition to developing simulation-based control capability, we plan to use linearity of the underlying equations and the explicit solvability property to develop efficient solutions extending classic Linear-Quadratic-Gaussian (LQG)-control approaches accounting for fluctuations. We will also explore respective robust formulations accounting for uncertainty in the system parameters.

• 3.

Short- and Long- term planning, aware of operations. The main problem in planning reliably over a sufficiently long time horizon (year and beyond) is in accounting for fluctuations and uncertainty of the heat consumption/production and also accounting for evolution and uncertainty in the network parameters (e.g. wear and tear of pipes and other devices) year(s) ahead. Combined with seasonal weather variability, this uncertainty, growing with time calls for including multiple consumption/production scenarios into a multi-level optimization designed to solve the planning problems. Moreover, to resolve operational quality of service constraints, subject to delays and inertia, plausible scenarios need to account for dynamics on the scale of hours to days.

• 4.

Including Demand Response (DR) elements in the district heating modeling and control. In this manuscript we have modeled heat consumption according to the oversimplified linear relation between the consumed heat and the temperature drop across the device. New concepts of DR, originally developed for electric loads [32, 33] but obviously extendable to heat loads, will require more general and flexible modeling of the heat consumption, which may be nonlinear, nonlocal in time and also changing with time depending on both exogenous an endogenous (control) signals.

• 5.

Integrated modeling and control of energy infrastructures. The concept of smart grids, when extended to integrated energy systems at the scale of a district or a city, has a lot of potential. To implement the concept we will attempt to take advantage of relations, dependencies and inter-operability of the power, natural gas and district heating networks. See, e.g., a recent example of a conceptually new take at the inter-operability of energy infrastructures considered in [34], where a unifying view on photovoltaic/natural gas power and cooling systems was proposed. To model joint operation of the three major energy infrastructures one will need to build and integrate dynamic models for all the relevant variables over the range of time scale starting from minutes and extending to days (and possibly even beyond). The technique developed in this manuscript to model heat flows in the district heating system is expected to contribute dynamic modeling, control and planning of the future integrated energy systems.

Appendix A An alternative Boundary Condition at the heat injection inlet of the pipe

As mentioned briefly in the main text, the choice of the Neumann, i.e. constant temperature, BC at the outlet of a pipe (or a bit more accurately at the hypothetical infinite outlet) is unambitious. However, selection of the BC at the heat injection, inlet, side of the pipe reflecting the nature of the heat injection is much less obvious. We suggest, based on physics considerations, that choosing the mixed BC corresponding to zero heat flux at the inlet is appropriate for representing a perfect heat source. This is simply because a good (ideal) heat source is build with the purpose of injecting the prescribed/controled amount of heat. However to reach an univocal conclusion on if the heat source model represents reality we will need to verify it through comparison with actual real world experiment on an operational heat source. This type of experimental test is out of scope of the paper. However, to provide input in this important line of future research we have experimented (in simulations) with different type of BC on the inlet (Neumann, Dirichlet or what we believe is the correct one - heat-flux constant) and have concluded that even thought the change in the inlet BC leads to some variability in a vicinity of the inlet, also proliferating along the pipe through an overall rescaling, the main physical picture of the front propagation and spreading remains invariant with respect to the choice of the inlet BC.

Therefore, in this Appendix we discuss, for completeness, an alternative choice of the BC at the inlet (15), the Neumann BC

[TABLE]

and show that the change in the inlet BC preserves solvability of the single pipe dynamics, critical for scalability of the large system analysis.

Analyzing Eqs. (14,25,16,17) via Laplace transform (in time) one derives

[TABLE]

Consider the following two, ”linear”-profile and ”exponential”-profile exemplary models of $\Delta(t)\doteq T^{(in)}(t)$

[TABLE]

where $\tau>0$ and $\Delta_{s}$ is the respective Laplace transform.

It is convenient to seek for the general solution of Eqs. (26) in the form

[TABLE]

where thus $\tilde{T}_{s}(x)$ represents Laplace transform of the deviation of temperature from the initial stationary solution.

The resulting equation for $\tilde{T}_{s}(x)$ and respective BC become

[TABLE]

thus producing the following explicit solution

[TABLE]

Then the complete solution (in the space-time domain) is presented in quadratures via the inverse Laplace transform

[TABLE]

According to Eqs. (37,38), $\tilde{T}_{s}(x)$ is an analytic function of complex $s$ anywhere in the complex plain except of the singularities (typically poles) associated with $T_{s}^{(in)}$ and also the square root branch cut singularity along the $]-\infty,s_{*}-\alpha-V^{2}/(4D)]$ cut , where

[TABLE]

Shifting the integration contour in Eq. (38) accordingly, evaluating the pole integrals and transforming the contour integral one arrives at the following explicit expression for the solution in quadratures (i.e. stated as an single-variable integral) in the case of exponential transient

[TABLE]

The remaining integral is of a Gaussian (diffusive) kernel type, which may thus be useful to represent as a series. Expending $\Delta_{s_{*}-Dq^{2}}$ in the Taylor series over $(-Dq^{2})$

[TABLE]

and evaluating the resulting integrals in Eq. (41) expicitly, one arrives at

[TABLE]

where $H_{n}(y)$ is the Hermit polynomial and $\ {}_{1}F_{1}(a,b,z)$ is the Kummer confluent hypergeometric function.

In the case of linear transient one derives

[TABLE]