**1. Introduction**

In the next few decades, the district heating (DH) will undergo substantial changes. The major operational changes associated with the 4th generation of DH are lower temperatures in distribution grids (the heat-carrier will be liquid water), enhancement of accumulation, and use of heat pumps or other temperature boosters [1,2]. Utilization of renewable sources and waste heat (both often fluctuate) aim at the reduction of carbon dioxide emissions and fossil fuel consumption [3]. Consumers, suppliers, and network providers will become members of the multilevel heat market. The production and consumption of heat or cold can be shifted between seasons using large long-term thermal storage, which has a lower relative price than smaller ones [4]. In such a complicated arrangement, the network must be sufficiently flexible, predictive, and intelligent.

Advanced and holistic operational optimization of the control strategies will become crucial. Although the operational optimization is usually being performed on already existing grids with given topologies, the simultaneous optimization of both grid topology and operation may be more beneficial than a separate approach. In other words, some grid topologies can allow better control strategies while others do not have to. The ideal control system should take long-term consequences into account (delays are most often affected by thermal inertia) [5]. The controller has to find the control actions by minimizing the sum of all the penalties (energy losses, demand mismatches, temporary price, etc.) defined inside the models of the individual components. This is called operational optimizations. It is crucial to fit as many simulations as possible to the control horizon (several minutes) to find the optimal values of the discrete control variables for the prediction horizon (several hours, days, weeks or even months). This depends on the models being fast and accurate. Regarding the pipes themselves, the correct prediction of step change propagation is most important as it determines the delay in temperature delivery. These delays are not only caused by pure advection through the system.

To ensure the sufficient performance of the model, the relative simulation time (ratio of physical time and computational time) is to be reduced. This might be achieved by a combination of several principles, such as lowering the number of used spatial derivatives (e.g., reducing to the one-dimensional description) or avoiding some of the discretizations altogether by using the analytic description of heat transfer [6] and friction [7] for all regimes in the pipes. When it comes to the hydraulic calculations, which are necessary to evaluate mass flows of heat carrier in the branches, it is usually considered satisfactory to use the static model, since the pressure waves spread significantly faster than thermal waves. It is also helpful to derive specific solutions that take advantage of known facts rather than those that cover a wide range of possible options (e.g., the presumption of incompressible fluid simplifies the dealing with the mass conservation laws significantly).

In previous work, a few main approaches that deal with simulations of DH grids are utilized. The first group is based on the presumption that the grid operates mostly in steady-state. Advanced techniques based on this can be found at [8], where authors develop a model and a method for parameter calibration that allows to better fit the measured data to simulation.

Another grea<sup>t</sup> example of steady-state estimation of DH grid is [9]. In this paper, the authors use automated customer meter system to evaluate mass flows, temperatures and consequently heat losses in pipes of a DH network. The optimization algorithm is utilized to estimate temperatures in the pipes, so simulated temperatures in the customer's location fit the measured data.

The second group of models, which include dynamic effects to some extent, is based on highly simplified models that directly link inlet temperature signal to the outlet signal. In the so-called node method [10], nodes that are connected by pipes of known geometrical and thermal properties represent the network. Based on mass flows in the pipes and the temperatures of water in individual nodes, the method steps through the time to evaluate temperature in all the nodes for all the time steps. It utilizes known time delays caused by advection and thermal inertias of the pipes.

Authors of [11] develop new so-called function method for the evaluation of thermal transients in DH grids. The method is based on the insertion of a Fourier series into the simplified Partial Differential Equation (PDE) of the fluid and solid regions of the circular pipe to derive the relationship between inlet and outlet temperature. Time delay and relative attenuations are incorporated into the model. The new method is reported to be about 37% faster than the standard node method. Model is validated in real DH system. No axial turbulent diffusion is considered.

Pure advection of inlet time of the water to the pipe is used in [12]. This value is then subtracted from the time when the liquid leaves the pipe. The resultant value represents the time the liquid spends inside the pipe which is then used for estimating the temperature drop due to heat loss. To incorporate thermal inertia, equivalent mixing volume is placed at the outlet. This results in a so-called plug-flow model. The model accurately simulated 168 h of operation of district heating grid in Pongau (Austria) in less than 2 s. No axial diffusion is incorporated.

Authors of [13] focus on transients in pipes with stable mass flows. They use steady-state axial temperature distributions and convolution with the axial diffusion kernel to arrive at an explicit mathematical function that describes the temporary transient behavior. Their solution is then compared with the numerical solution of the appropriate Partial Differential Equation (PDE) with positive results. No dynamic re-accumulation in the pipe walls is included.

A steady-state heat transfer model is combined with a variable transport delay model in [14] to obtain a fast model capable of simulations with variable flows. The model resulted in approximately 4000 times less intensive computation in comparison with one-dimensional PDE approach (both models were developed in Matlab/Simulink). This speedup factor rises with the utilization of larger time steps but at the cost of higher error. Axial diffusion and thermal inertia of the wall is considered.

Comparison of the node method and pseudo transient approach with experimental data is conducted in [15]. Authors of the article conclude that thermal diffusion caused by turbulent behavior of the fluid inside the pipe has a significant effect on the thermal wave propagation. The effect of turbulent diffusion is more significant with sharper inlet temperature changes and that the proper modeling should consider this. Further, they conclude that the pseudo transient approach predicts the thermal propagation more accurately than the node method.

The third group of models is based on the discretization of the PDE describing the pipe in one or more dimensions (numerical/physical models). Method of characteristics with the third-order numerical scheme is used in [16] to solve the energy Equation. It seems to preserve the shapes and amplitudes of the heat wave, which indicates that there is little or none numerical diffusion in the solution. Axial diffusion and thermal dynamics of the wall is not incorporated.

Improved third order numerical solution with total variation diminishing property (reduces oscillations near sudden temperature changes) is used in [17]. Simplified dynamic accumulation of heat into the pipe's wall is considered. An indicator that describes the influence of the thermal inertia of the wall onto the wave propagation is developed. The proposed model is validated in a real system. Axial turbulent diffusion is not incorporated into the model.

Although numerical models are usually regarded as slow, they have one significant advantage. They can capture most of the dynamic phenomena affecting the traveling thermal wave in considerable detail. According to the review of relevant literature given above, any of the states of the art models (of any group) does not incorporate all the mentioned effects at once. To speed up numerical models, it is desirable to limit the number of operations as well as the number of memory allocations (saving fewer data and creating fewer temporaries). To obtain a smooth solution, there is a need to choose minor steps unless an appropriate nonlinear interpolation is possible.

Usually, problems, such as advection, are numerically solved by finite volume methods (FVM). Either explicit or implicit schemes are used. The PDE is often discretized along all dimensions, including the time. The quadratic upstream interpolation with estimation terms (QUICKEST) is often considered to be the most efficient scheme [18–21]. A good FVM scheme should preserve sharp temperature changes while not introducing unphysical oscillation. Non-oscillatory schemes are nonlinear and usually involve a flux limiter, which makes them more expensive during run-time. The only linear non-oscillatory FVM is first order upwind scheme, but it suffers from severe unphysical (numerical) diffusion. Luckily, when real diffusion is involved, the oscillation of higher-order schemes fades out. An example of this is presented in the following section.

Heat transfer is very often solved by finite element method (FEM), which is based on a weak formulation of the PDEs. Nodal values of the stencil amplify local non-zero normalized field functions (shape functions). The global temperature field function is then a superposition of the magnified local fields. Therefore, a good FEM scheme is produced by shape functions that closely resemble the true local temperature distribution. In the following section, there is a detailed description of an approach that results in the stencil coefficients for the radial model of a solid wall. Steady-state biased shape functions are used to capture the less transient periods very accurately even through roughest meshes.

In order to facilitate the various models of individual components of a DH system, all components should be solved simultaneously to include their mutual effects. Method of lines can be utilized here. In this method, the PDEs are approximated by sets of ordinary differential equations (ODE), where only the spatial derivatives are discretized. Such models are easy to combine using an ODE solver.

Open-source environments with powerful ODE solvers can be used to solve these systems. Two interesting ones are the OpenModelica and Julia languages. The first uses a component-oriented, Equation-based language called Modelica. It specializes in the modeling of complex systems dealing with various physics backgrounds [22]. There are many libraries containing predefined models of components which can be easily combined together to create a complex system. The OpenModelica

compiler compiles code that solves efficiently time-dependent variables based on given parameters. It uses solvers for Differential algebraic Equations (DAE) to solve implicit systems, which allows the user to specify the relationship between variables without the need for being explicit. Communication with MATLAB can be also established through scripting and file exchange [23] (this is advantageous because MATLAB's default ODE solvers are very slow).

Julia is a general-purpose, high-level, dynamic language that approaches the speeds of statically-compiled languages like C [24]. There are several tools for benchmarking and profiling of written algorithms, which helps to identify their potential speedups. There are many well-developed packages already, amongs<sup>t</sup> which is a very performant (meaning that it is working well enough to be considered functional but may not ye<sup>t</sup> be fully optimized) package dealing with Differential Equations [25]. This package contains many solvers and can even interface to advanced C libraries, such as Sundials [26]. The results use interpolation with predefined tolerance to return values at an arbitrary time. This together makes the language highly suitable for scientific computing.

The aim of this paper was to derive and test specialized and fast numerical model of a simple DH pipe that incorporates sharp interface preservation, turbulent diffusion, and detailed radial thermal inertia evaluation caused by the wall.

#### **2. Materials and Methods**

This chapter contains a very in-depth description of the mathematics behind the model. In order to ensure that the resultant model is fast, a one-dimensional PDE is used to avoid two spatial derivatives. Here, the incompressibility of the fluid is assumed. Individual variables depend on the position and time as denoted by variables *x* and *t* in the following equation:

$$\frac{dT(\mathbf{x},t)}{dt} = -u(t)\frac{dT(\mathbf{x},t)}{d\mathbf{x}} + D\frac{d^2T(\mathbf{x},t)}{d\mathbf{x}^2} + \mathcal{S}(\mathbf{x},t),\tag{1}$$

where *T* is the local cross-sectional temperature of the fluid, *u* is the axial fluid velocity, *D* is the axial diffusion coefficient, and *S* is a source term (dealing with heat flow from inner wall and friction). Because the method of lines is used here, only the spatial derivatives must be discretized. The pipe is divided into *na* cylindrical volumes, where spatially dependent variables are localized (see Figure 1). Localized temperatures represent the volume-averaged temperature and, localized sources represent the incoming heat from the solid shell divided by the heat capacity of one volume element.

**Figure 1.** Division of pipe's fluid region into volumes, with localized spatial dependent variables.

Finite differences are used for discretization. There is a simple relationship between the finite difference schemes for any derivation order (including zero) and given *n* stencil points [27]:

$$
\varepsilon = \frac{1}{\Delta x^{\text{m}}} Z^{-1} \delta \,, \tag{2}
$$

where *c* is a vector of the stencil coefficients, *m* is the derivation order, *Z* is an *n* × *n* matrix such as:

$$Z\_{i,j} = \left(\frac{\mathbf{x}\_i - \mathbf{x}\_0}{\Delta \mathbf{x}}\right)^{j-1},\tag{3}$$

in which *xi* is the axial position of a stencil point, *x0* is the axial position of the point in which the derivative is demanded, and *δ* represents a vector such as:

$$\delta\_{\vec{i}} = \begin{cases} 0 & \text{ $\vec{i} \neq m+1$ } \\ m! & \text{ $\vec{i} = m+1$ } \end{cases} \tag{4}$$

The first term on the right-hand side of Equation (1) dealing with advection is discretized in two ways. The first way utilizes Equation (2) directly, which results in the finite difference method (FDM). The second way is based on FVM, where Equation (2) is used for the evaluation of the values at element faces (zero order derivatives). The discretization of the advection part should be good in shock capturing; therefore, 660 tests are conducted to compare this ability among the different schemes. By filling rows with the stencil coefficients, a script written in Julia constructs the matrix *A* that approximates the derivative operator. Among the input parameters, there is a number of used upwind and downwind stencil points. While filling the first few rows of matrix *A,* both these parameters are automatically decreased. When filling the last few rows, only the downwind parameter is decreased because of the missing nodes. For example, the maximum number of available upwind nodes at the inlet is one (the boundary temperature). No oscillation is caused by central schemes near the boundaries (unless a central scheme is used explicitly) because the local upwind bias is always preserved (first rows) or increased (last rows). The boundary vector *Ain* is constructed as well. Both are used to evaluate the pure advection problem in the following form:

$$\frac{dT}{dt} = -u(AT + T\_{\text{in}}A\_{\text{in}})\_{\text{\textdegree}} \tag{5}$$

where *T* is the vector of localized fluid volume temperatures and *Tin* is the time-dependent scalar representing the inlet (boundary) temperature signal. Problems are solved using normalized variables, and the behavior of the schemes is presented in a normalized variable diagram (NVD). Sundials library is used for solving the problem in time, namely the CVODE algorithm with backward differencing and with relative and absolute tolerances of 10−6. It has an adaptive time step and produces a continuous temperature signal at the outlet by using Hermite interpolation. Examples of NVDs are shown in Figure 2. The diagrams contain hints about the schemes (*U* and *D* are the number of upwind and downwind nodes, respectively, and *n* is a number of volume elements).

**Figure 2.** Examples of bad-behaving advection schemes. (**a**): Low order upwind biased scheme causes a lot of numerical diffusion; (**b**): Low order central scheme is being too oscillatory.

It is observed that central schemes have a stronger tendency to oscillate (Figure 2b) and are more computationally expensive than upwind-biased schemes (see Figure 3). Overall, higher order upwind schemes (see Figure 4) seem to show a good compromise between the degree of their oscillation, their ability for shock capturing, and computational expensiveness.

**Figure 3.** Comparison of computational time of the tests between upwind-biased and central schemes for FDM. Similar results can be observed for FVM. The values in brackets represent rounded row stencil coefficients in the typical row of the derivative operator *A* (when Δ*x* = 1).

**Figure 4.** Examples of well-behaving advection schemes. (**a**): Higher order upwind biased FVM scheme with mild numerical diffusion and mild oscillations; (**b**): Higher order upwind biased FDM scheme with mild numerical diffusion and very mild oscillations.

The first order central scheme is used for the construction of *B*, which is the matrix that approximates the second derivative operator dealing with axial diffusion in Equation (1). Boundary vector *Bin* allows the calculation of the second derivative at the inlet node as well. Equation (5) is extended and results in Equation (6):

$$\frac{dT}{dt} = -\mu(AT + T\_{\text{in}}A\_{\text{in}}) + D(BT + T\_{\text{in}}B\_{\text{in}}),\tag{6}$$

diffusion coefficient *D* may either be lumped or localized into the form:

$$D\_i = w\_{1,i} \cdot u \cdot d\_1 + w\_{2,i} \tag{7}$$

The first term in Equation (7) deals with turbulent diffusion. The second term, which is strictly non-negative, corresponds to velocity-independent effects, such as conductive diffusion (which is shown later to be negligible) or spreading due to effects of gravity. If the second term in Equation (7) is neglected, the Equation (6) can be simplified:

$$\frac{dT}{dt} = u\left[ (d\_1 B - A)T + (d\_1 B\_{in} - A\_{in})T\_{in} \right] = u(MT + M\_{in}T\_{in}),\tag{8}$$

where *M* and *Min* are the constant matrix operator and boundary vector of the fluid region, respectively. This step should improve the performance of the model. It also makes the Peclét number independent of flow velocity. This means that when proper axial discretization is chosen to ensure non-oscillatory behavior of the method, it is preserved for all velocities of the fluid. Automatic axial discretization, based on values *w1,i*, is then possible. The weights in Equation (7) allow the tuning of the pipe to mimic a real installation more accurately (similar weighting is possible for the source term as well).

Using the analytic solution to step input in the NVD, the effect of *Di* on oscillation of Equation (6) is studied. Any value of the numerical solution that is outside of the unit interval is considered to be an overshot caused by oscillation. This is a very simple way to estimate the residual oscillation of Equation (6). FDM is chosen for the model because of its lower oscillations (see Figure 5).

**Figure 5.** Effect of physical diffusion on unphysical oscillation on sharp interfaces for *na* = 100 (the critical Peclét number, which makes the residual oscillation less than 10−6. It is approximately 6.47 for the U2D1FVM and 8.51 for the U2D1FDM).

The last term of Equation (1) is more complicated since it must deal with heat dynamics in the solid shells surrounding the pipe. From the fluid viewpoint, the only important value is the inner wall temperature of the first shell. Using publication [6], the Nusselt number for all flow regimes is:

$$Nu = \begin{cases} 3.66 \text{ if } Re \le 2300 \\ 3.52 \cdot \mathcal{U}^4 - 45.148 \cdot \mathcal{U}^3 + 212.13 \cdot \mathcal{U}^2 - 427.45 \cdot \mathcal{U}^2 + 316.08 \text{ if } 2300 \le Re \le 3100 \\\ \frac{f\_r}{8} \cdot \frac{(Re - 1000) \cdot \mathcal{Re}}{1 + 12.7 \cdot \left(\frac{f\_r}{8}\right)^{1/2} \cdot \left(Pr^{2/3} - 1\right)} \text{ if } Re < 3000 \end{cases},\tag{9}$$

where *Re* is the Reynolds number, *U* = *Re/1000*, *fr* is the Darcy-Weisbach friction factor, and *Pr* is the Prandtl Number. The friction factor might be evaluated using the Churchill Equation presented in [7] for all flow regimes:

$$f = 8 \cdot \left( \left( \frac{8}{Re} \right)^{12} + \frac{1}{\left( \theta\_1 + \theta\_2 \right)^{3/2}} \right)^{1/12},\tag{10}$$

where

$$\begin{array}{c} \theta\_1 = \left(-2.457 \cdot \log\left(\left(\frac{7}{R\varepsilon}\right)^{0.9} + 0.27 \cdot \left(\frac{\varepsilon}{d\_1}\right)\right)\right)^{16} \\ \theta\_2 = \left(\frac{37530}{R\varepsilon}\right)^{16} \end{array} \tag{11}$$

in which *ε* is the roughness of the inner wall and *d1* is the inner diameter. It is a simple step to evaluate the local source term *S* in the form of a vector:

$$S\_i = \frac{4 \cdot Nu \cdot \lambda}{d\_1^2 \cdot \rho \cdot c\_p} \cdot (T\_{sn,i,1} - T\_i). \tag{12}$$

where *Tsn,i,1* represents the inner wall (later first node of the solid region) temperatures of the solid shell, *cp* is the specific heat capacity of the fluid, *ρ* is the mass density of the fluid, and *λ* is the heat conductivity of the fluid.

The time-dependent distribution of temperature in the solid shells affects the axial transients of the heat waves (charging and discharging of the accumulation layers). Tangential and axial heat conductions in the shells are ignored because even when the flow velocity is small, the rate in which the heat is being transported by advection is often much higher than that by conduction. This might be expressed in the form of the Péclet number for conduction in both the fluid:

$$Pe\_f = d\_1 \cdot u \cdot \frac{\rho \cdot \varepsilon\_p}{\lambda} \approx d\_1 \cdot u \cdot 6,\\ 29 \cdot 10^6 \gg 1,\tag{13}$$

and the solid (e.g., steel):

$$Pe\_s = d\_1 \cdot u \cdot \frac{\rho\_s \cdot \mathcal{c}\_{p,s}}{\lambda\_s} \approx d\_1 \cdot u \cdot 8 \cdot 10^4 \gg 1. \tag{14}$$

As shown, the inner diameters of the pipes or the flow velocities would have to be very small for the conductive axial diffusion to be comparable with advection. This is not common with real applications in the DH systems which also reduces the number of necessary spatial derivatives used in the model. The fluid volumes can be mathematically coupled with one-dimensional models of purely radial heat transfer (later solid segments) described by PDE in the form:

$$\frac{dT\_{s,i}}{dt} \cdot \rho\_s \cdot c\_{ps} = \frac{1}{r} \frac{d}{dr} \left(\lambda\_s \cdot r \cdot \frac{dT\_{s,i}}{dr}\right),\tag{15}$$

where *r* is the radial coordinate, *Ts,i* is the continuous temperature in the solid segment, *ρs* is the mass density, *cps* is the specific heat capacity, and *λs* is the heat conductivity. Equation (15) might be rewritten into the weak form using the shape function *V* (notice that integrated term in the bracket affects the system only at the boundaries).

$$2\cdot \pi \cdot \Delta \mathbf{x} \cdot \int\_{r\_1}^{r\_2} \mathbf{r} \cdot V \cdot \frac{d\mathbf{T}\_{sj}}{dt} \cdot \rho\_\delta \cdot \mathbf{c}\_{\mathrm{pS}} \cdot dr = 2 \cdot \pi \cdot \Delta \mathbf{x} \cdot \left[ -\int\_{r\_1}^{r\_2} \mathbf{r} \cdot \frac{dV}{dr} \cdot \frac{dT\_{sj}}{dr} \cdot \lambda\_\delta \cdot dr + \left[ \mathbf{r} \cdot V \cdot \lambda\_\delta \cdot \frac{dT\_{sj}}{dr} \right]\_{r\_1}^{r\_2} \right]. \tag{16}$$

The steady state temperature distribution is used as the base for the nodal shape functions:

$$T\_{s\varepsilon,i,j}(r) = \frac{\ln(r/\mathcal{R}\_{\hat{\gamma}}) \cdot \left(T\_{\alpha i,j+1} - T\_{\alpha i,j}\right)}{\ln\left(R\_{\hat{\gamma}+1}/R\_{\hat{\gamma}}\right)} + T\_{s\alpha j,j} = T\_{s\alpha j,j+1} \cdot \frac{\ln\left(r/\mathcal{R}\_{\hat{\gamma}}\right)}{\ln\left(R\_{\hat{\gamma}+1}/R\_{\hat{\gamma}}\right)} + T\_{s\alpha j,j} \left(1 - \frac{\ln\left(r/R\_{\hat{\gamma}}\right)}{\ln\left(R\_{\hat{\gamma}+1}/R\_{\hat{\gamma}}\right)}\right),\tag{17}$$

where *Tse,i,j* are the temperature distributions in the cylindrical elements, *Tsn,i,j* are nodal temperatures, and *Rj* are radial coordinates of the nodes. Equation (17) is the solution to the following PDE with given boundary conditions (see also Figure 6):

$$\begin{array}{c} \frac{1}{r} \frac{d}{dr} \left( \lambda\_{s,j} \cdot r \cdot \frac{dT\_{sn,j}}{dr} \right) = 0\\ T\_{sc,i,j}(R\_{j+1}) = T\_{sn,i,j+1} \quad \text{ } \tag{18} \\ T\_{sc,i,j}(R\_j) = T\_{sn,i,j} \end{array} \tag{18}$$

**Figure 6.** Graphical representation of the variables.

Nodal shape functions are derived from Equation (17) by considering the effect of one node to the global temperature distribution. Therefore, they take the following piecewise form (notice that boundary nodes are neighboring one element only, so the shape functions have only the appropriate half):

$$V\_{sn,j}(r) = \begin{cases} \begin{cases} & 1 - \frac{\ln\left(r/R\_j\right)}{\ln\left(R\_{j+1}/R\_j\right)}, R\_j \le r \le R\_{j+1} \\ & 0, R\_j > r > R\_{j+1} \\ & \frac{\ln\left(r/R\_{j-1}\right)}{\ln\left(R\_j/R\_{j-1}\right)}, R\_{j-1} \le r < R\_j \\ & 1 - \frac{\ln\left(r/R\_j\right)}{\ln\left(R\_{j+1}/R\_j\right)}, R\_j \le r \le R\_{j+1} \\ & 0, R\_{j-1} > r > R\_{j+1} \\ & \frac{\ln\left(r/R\_{j-1}\right)}{\ln\left(R\_j/R\_{j-1}\right)}, R\_{j-1} \le r < R\_j \\ & 0, R\_{j-1} > r > R\_j \end{cases}, j = n$$

When mass densities and specific heat capacities of the elements are constant, the following integration (according to the left-hand side of Equation (16)) results in discretized heat capacity:

$$\mathbf{C}\_{\rm sn,j} = 2 \cdot \boldsymbol{\pi} \cdot \Delta \mathbf{x} \cdot \left[ \rho\_{\rm s,j-1} \cdot \mathbf{c}\_{\rm ps,j-1} \cdot \int\_{R\_{j-1}}^{R\_j} \boldsymbol{r} \cdot \mathbf{V}\_{\rm sn,j} \cdot d\mathbf{r} + \rho\_{\rm s,j} \cdot \mathbf{c}\_{\rm ps,j} \cdot \int\_{R\_j}^{R\_{j+1}} \boldsymbol{r} \cdot \mathbf{V}\_{\rm sn,j} \cdot d\mathbf{r} \right],\tag{20}$$

where *Csn,j* are the nodal heat capacities, *ρse,i,j* are the mass densities of the elements, and *cps,j* is the specific heat capacity of the elements. For inside nodes, the integration of (20) results in:

$$\mathbf{C}\_{\text{sti},\boldsymbol{\zeta}} = \boldsymbol{\pi} \cdot \Delta \mathbf{x} \cdot \left[ \rho\_{s\_{\boldsymbol{\zeta}} - 1} \cdot \mathbf{c}\_{p s\_{\boldsymbol{\zeta}} - 1} \cdot \left( \frac{\boldsymbol{R}\_{j-1}^2 - \boldsymbol{R}\_j^2}{2 \ln \left( \boldsymbol{R}\_{\boldsymbol{\zeta}} / \boldsymbol{R}\_{\boldsymbol{\zeta} - 1} \right)} + \boldsymbol{R}\_j^2 \right) - \rho\_{s\_{\boldsymbol{\zeta}}} \cdot \mathbf{c}\_{p s\_{\boldsymbol{\zeta}}} \cdot \left( \frac{\boldsymbol{R}\_j^2 - \boldsymbol{R}\_{j+1}^2}{2 \ln \left( \boldsymbol{R}\_{\boldsymbol{\zeta} + 1} / \boldsymbol{R}\_{\boldsymbol{\zeta}} \right)} + \boldsymbol{R}\_j^2 \right) \right],\tag{21}$$

*Energies* **2019**, *12*, 633

For the boundary near the fluid region, there is only the second term in the bracket (there is no element below), while for the outer boundary node, there is only the first term (there is no element above).

The net heat flow incoming to the nodes in the solid region is (according to the right-hand side of Equation (16)):

$$H\_{\rm sn,j,j} = 2 \cdot \pi \cdot \Delta \mathbf{x} \cdot \left[ T\_{\rm sn,j,j-1} \cdot h\_{\rm sn,j-1,j} + T\_{\rm sn,j,j} \cdot h\_{\rm sn,j,j} + T\_{\rm sn,j,j+1} \cdot h\_{\rm sn,j+1,j} \right],\tag{22}$$

where

$$h\_{\mathbf{s}n,j-1,j} = -\int\_{R\_{j-1}}^{R\_j} r \cdot \lambda\_{\mathbf{s},j-1} \cdot \frac{dV\_{\mathbf{s}n,j-1}}{dr} \frac{dV\_{\mathbf{s}n,j}}{dr} dr = \frac{\lambda\_{\mathbf{s},j-1}}{\ln(R\_j/R\_{j-1})}$$

$$h\_{\mathbf{s}\mathbf{s},j,j} = -\int\_{R\_{j-1}}^{R\_j} r \cdot \lambda\_{\mathbf{s},j-1} \cdot \left(\frac{dV\_{\mathbf{s}n,j}}{dr}\right)^2 dr - \int\_{R\_j}^{R\_{j+1}} r \cdot \lambda\_{\mathbf{s},j} \cdot \left(\frac{dV\_{\mathbf{s}n,j}}{dr}\right)^2 dr = \frac{-\lambda\_{\mathbf{s},j-1}}{\ln(R\_j/R\_{j-1})} + \frac{-\lambda\_{\mathbf{s},j}}{\ln(R\_{j+1}/R\_j)}.\tag{23}$$

$$h\_{\mathbf{s}n,j,j} = -\int\_{R\_j}^{R\_{j+1}} r \cdot \lambda\_{\mathbf{s},j} \cdot \frac{dV\_{\mathbf{s}n,j+1}}{dr} \frac{dV\_{\mathbf{s}n}}{dr} dr = \frac{\lambda\_{\mathbf{s},j}}{\ln(R\_{j+1}/R\_j)}$$

The nodal net heat flow changes for boundary nodes. For the inner boundary node (here the integrated term in brackets of Equation (16) is replaced by heat flow from fluid to solid), the net heat flow is:

$$H\_{\rm sn,j,1} = 2 \cdot \pi \cdot \Delta \mathbf{x} \cdot \left[ T\_{\rm sn,j,1} \cdot h\_{\rm sn,1,1} + T\_{\rm sn,j,2} \cdot h\_{\rm sn,2,1} \right] + \pi \cdot \Delta \mathbf{x} \cdot \mathbf{N} \boldsymbol{u} \cdot \boldsymbol{\lambda} \cdot \left( T\_{\rm i} - T\_{\rm sn,j,1} \right), \tag{24}$$

where

$$\begin{aligned} h\_{\text{sn},1,1} &= -\int\_{R\_1}^{R\_2} r \cdot \lambda\_{s,1} \cdot \left(\frac{dV\_{\text{sn},1}}{dr}\right)^2 dr = \frac{-\lambda\_{s,1}}{\ln(R\_2/R\_1)}\\ h\_{\text{sn},2,1} &= -\int\_{R\_1}^{R\_2} r \cdot \lambda\_{s,1} \cdot \frac{dV\_{\text{sn},2}}{dr} \frac{dV\_{\text{sn},1}}{dr} dr = \frac{\lambda\_{s,1}}{\ln(R\_2/R\_1)} \end{aligned} \tag{25}$$

and similarly for the outer boundary node:

$$H\_{\rm sn,i,n} = 2 \cdot \pi \cdot \Delta \mathbf{x} \cdot \left[ T\_{\rm sn,i,n-1} \cdot h\_{\rm sn,n-1,n} + T\_{\rm sn,i,n} \cdot h\_{\rm sn,n,n} \right] + a\_{\rm out} \cdot 2 \cdot \pi \cdot R\_{\rm n} \cdot \Delta \mathbf{x} \cdot \left( T\_{\rm air} - T\_{\rm sn,i,n} \right), \tag{26}$$

where the last term represents true heat loss (heat transfer to the environment), *αout* is the convective heat transfer coefficient at the outer wall, *Tair* is the ambient temperature and:

$$\begin{split} h\_{\text{sII},\text{n}-1,\text{n}} &= -\int\_{R\_{n-1}}^{R\_{\text{n}}} r \cdot \lambda\_{s,\text{n}-1} \cdot \frac{dV\_{\text{s},\text{n}-1}}{dr} \frac{dV\_{\text{s},\text{n}}}{dr} dr = \frac{\lambda\_{s,\text{n}-1}}{\ln(R\_{\text{n}}/R\_{\text{n}-1})} \\ h\_{\text{sII},\text{n},\text{n}} &= -\int\_{R\_{\text{n}-1}}^{R\_{\text{n}}} r \cdot \lambda\_{s,\text{n}-1} \cdot \left(\frac{dV\_{\text{s},\text{n}}}{dr}\right)^{2} dr = \frac{-\lambda\_{s,\text{n}-1}}{\ln(R\_{\text{n}}/R\_{\text{n}-1})} \end{split} \tag{27}$$

Using Equations (16), (21) and (22), the set of ODEs describing heat dynamics in the solid is:

$$\frac{dT\_{sn,i,j}}{dt} = \frac{H\_{sn,i,j}}{C\_{sn,j}},\tag{28}$$

Values in matrix *hsn,i,j* are nonzero only for values given by Equations (23), (25) and (27). Under such conditions, the following matrices might be constructed:

$$k\_{i,j} = 2 \cdot \pi \cdot \Delta x \cdot \frac{h\_{sn,i,j}}{\mathbb{C}\_{sn,j}} \, ^\prime \tag{29}$$

*Energies* **2019**, *12*, 633

$$k\_{\rm in,i,j} = \begin{cases} \displaystyle \;j = 1, \; \pi \cdot \Delta \mathbf{x} \cdot \mathbf{N} \mathbf{u} \cdot \boldsymbol{\lambda} \cdot \left(T\_{\bar{i}} - T\_{\rm sn,i,1}\right) / C\_{\rm sn,1} \\ \displaystyle \; \mathbf{2} \le j \le n - 1, \; \boldsymbol{0} \\ \displaystyle \;j = n, \; \boldsymbol{\alpha}\_{\rm out} \cdot \mathbf{2} \cdot \pi \cdot \mathcal{R}\_{\rm n} \cdot \Delta \mathbf{x} \cdot \left(T\_{\rm air} - T\_{\rm sn,i,n}\right) / C\_{\rm sn,n} \end{cases} \tag{30}$$

Using Equations (29) and (30), the relation (28) might be expressed in a form similar to Equation (6):

$$\frac{dT\_{sn,i}}{dt} = k \cdot T\_{sn,i} + k\_{in,i} \,. \tag{31}$$

where *Tsn,i* are vectors of the nodal temperatures of the solid segments, *k* is a constant matrix (constructed by Equation (29)) capturing the heat dynamics of the nodes inside segments, and *kin,i* are the time-dependent vectors containing boundary conditions of the segments. Banded matrices can be used in order to express the dynamics fully in a single equation:

$$\frac{dT\_{\rm sn}}{dt} = K \cdot T\_{\rm sn} + K\_{\rm in}.\tag{32}$$

where *Tsn* is the column vector filled with column vectors *Tsn,i*, *K* is the banded square matrix filled with *k* on its diagonal, and *Kin*is the time-dependent boundary vector filed with column vectors *ksn,i*.

This description of the solid region will satisfy the analytic solution completely at steady states but will have a tendency to deviate more with rapid changes of fluid temperature. Individual layers can be further divided into sublayers, which will increase the accuracy during nonstationary events.

#### **3. Results and Discussion**

The new model of radial heat transfer is tested for two different cases (see Table 1). The first case is not an actual installation but serves as an example of a case with very strong accumulation. A sudden change in the fluid temperature causes high temporal heat flux between the fluid and the wall. A quick reduction in the heat flux follows due to a rise of the innermost wall surface temperature. The radial elements increasing in size from the inner wall outward should capture this better. Smaller time-steps must be then used to ensure numerical stability. Element sizes are defined by a grow factor. Qualities of the solutions at the innermost piece of the wall are studied using a step change in fluid temperature from 0 to 10 K.


**Table 1.** The property material of the steel circular wall (inner pipe diameter is 0.6 m).

Since the temperature of the fluid steps up to 10 K and the correct inner wall temperature after 10 s is about 2.16 K, the error in heat flux caused by the use of 5 subdivisions (with grow factor 1) would be over 10%. Grow factor 1.0 with 10 subdivisions provides almost equivalent results to an advanced solution that uses grow factor 2.0 with 5 subdivisions, but it is slightly faster (see Figure 7). The thin wall in the second case behaves much better. It produces solutions accurate enough even with lower subdivisions (compare Figures 8–11 with Figures 12–15). This ensures good computational speed for such arrangements. To further estimate the accuracy of the new radial model for the case with strong accumulation, the maximal errors are estimated using temperature refinements between two consecutive levels of subdivision (see Figures 8–11). The whole event of heat recharging is used (7200 s). A tolerance of 10−<sup>12</sup> is chosen to ensure that the time stepping method can always take advantage of a superior spatial discretization.

**Figure 7.** Demonstration of the solution quality (in the radial model) at 10 s for different subdivisions (Case 1).

**Figure 8.** Absolute maximum temperature error estimations of the inner wall surface (Case 1).

**Figure 9.** Absolute maximum temperature error estimations of the outer surface (Case 1).

**Figure 10.** Relative maximum temperature error estimations of the inner surface (Case 1).

**Figure 11.** Relative maximum temperature error estimations of the outer surface (Case 1).

**Figure 12.** Absolute maximum temperature error estimations of the inner surface (Case 2).

**Figure 13.** Absolute maximum temperature error estimations of the outer surface (Case 2).

**Figure 14.** Relative maximum temperature error estimations of the inner surface (Case 2).

**Figure 15.** Relative maximum temperature error estimations of the outer surface (Case 2).

The errors in the second case for lower subdivisions are equivalent to errors in the first case for high subdivisions. This is because the inertia of the thin wall is not significant. The relative errors of the inner and outer walls are almost the same (Figures 13 and 14).

The composite model of the pipe written in OpenModelica is tested in a similar way (Julia was significantly slower here, possibly due to the use of numerical jacobians instead of symbolic ones). The 100 m long pipe with a cross-section corresponding to the first case of the previous test set is initialized with a temperature of 323.15 K. After steady-state is reached, a 30 K rise and drop in the inlet temperature generate the transient process. The duration of the tests was 16,000 s. The tolerance is set to <sup>10</sup>−6, and the values are saved every 3 s. Error estimations are evaluated using a discrete gradient (backward differencing) between the set of discretization. Using contours, the "willingness to refine maximal errors" of the model can be mapped (see Figures 16 and 17). The computational time can be mapped directly (see Figure 18).

**Figure 16.** Estimations of maximal absolute errors at the outlet of the pipe (Case 1), where the red arrow signifies the ideal direction for discretization.

**Figure 17.** Estimations of maximal relative errors at the outlet of the pipe (Case 1), where the red arrow signifies the ideal direction for discretization.

Similar tests for the second case clearly demonstrate the independence on the radial number of subdivisions (see Figures 19 and 20). There is nothing to be refined in the radial direction. Further subdivisions only take more computational time (see Figure 21). The higher computational time in the right top corner of Figure 21 is caused by smaller time steps to ensure numerical stability.

**Figure 18.** Computation times of the composite model tests (Case 1).

**Figure 19.** Estimations of maximal absolute errors at the outlet of the pipe (Case 2), where the red arrow signifies the ideal direction for discretization.

**Figure 20.** Estimations of maximal relative errors at the outlet of the pipe (Case 2), where the red arrow signifies the ideal direction for discretization.

**Figure 21.** Computation times of the composite model tests (Case 2).

The Figure 22 demonstrates the effect of the strong accumulation using the outlet temperature perspective. The above results deal with the model alone. Unfortunately, real data were not available to the authors to validate the model properly. The model is at least compared (partially validated) to a very fine model developed in Fluent (axisymmetric 2D in Figures 23 and 24) and STAR-CCM+ (fully 3D in Figures 25 and 26). The tolerance of the new model was lowered to 10−<sup>5</sup> and samples were saved every 10 s to demonstrate the effectiveness of the new model. The time span of the simulation was 16,000 s.

**Figure 22.** Comparison of Case 1 and Case 2 (with constant mass flow of 282.74 kg/s).

**Figure 23.** Comparison with a model from Fluent (the solid material property differ slightly here—mass density is 7850 kg/m3, heat conductivity is 50 W/(mK), and specific heat capacity is 500 J/(kgK)).

**Figure 24.** Detail of the difference in temperature on comparison with a model from Fluent.

**Figure 25.** Comparison with a model from STAR-CCM+.

**Figure 26.** Detail of the difference in temperature on comparison with a model from STAR-CCM+.

Although the maximal temperature difference is higher than the model from STAR-CCM+, the temperature difference falls close to zero quicker than the model from Fluent. The effect of gravity on wave propagation is such that it behaves as additional axial diffusion, which can be tuned up by adjusting weights in Equation (7). Figure 26 shows this fact by the more time-spread rise in temperature at the outlet. Figures 27 and 28 show how the gravity spreads the temperature front in the STAR-CCM+ model 7 s after both hot and cold waves entered the pipe.

**Figure 27.** Effect of gravity on the temperature distribution of a hot wave (the upper part shows the temperature front without effects of gravity, the lower part shows the contrary).

**Figure 28.** Effect of gravity on the temperature distribution of cold wave (the upper picture shows the temperature front without effects of gravity, the lower picture shows the contrary).

The comparison of the new model with the very fine numerical counterparts shows that the new model overestimates the wave to reach the outlet slightly later (asymmetric peaks in temperature difference in the Figures 23–26). This might be caused by the fact that the temperature and velocity profiles in the new model are considered to be fully developed right from the beginning, while with the very fine counterparts they are not.

#### **4. Summary of the Results**

