Performant and Simple Numerical Modeling of District Heating Pipes with Heat Accumulation

Abstract

This paper compares approaches for accurate numerical modeling of transients in the pipe element of district heating systems. The distribution grid itself affects the heat flow dynamics of a district heating network, which subsequently governs the heat delays and entire efficiency of the distribution. For an efficient control of the network, a control system must be able to predict how “temperature waves” move through the network. This prediction must be sufficiently accurate for real-time computations of operational parameters. Future control systems may also benefit from the accumulation capabilities of pipes. In this article, the key physical phenomena affecting the transients in pipes were identified, and an efficient numerical model of aboveground district heating pipe with heat accumulation was developed. The model used analytical methods for the evaluation of source terms. Physics of heat transfer in the pipe shells was captured by one-dimensional finite element method that is based on the steady-state solution. Simple advection scheme was used for discretization of the fluid region. Method of lines and time integration was used for marching. The complexity of simulated physical phenomena was highly flexible and allowed to trade accuracy for computational time. In comparison with the very finely discretized model, highly comparable transients were obtained even for the thick accumulation wall.

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 great 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,19,20,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, amongst which is a very performant (meaning that it is working well enough to be considered functional but may not yet 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\left(x,t\right)}{dt}=-u\left(t\right)\frac{dT\left(x,t\right)}{dx}+D\frac{d^{2}T\left(x,t\right)}{d{x}^2}+S\left(x,t\right),$$

(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 n_a 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.

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]:

$$c=\frac{1}{\Delta x^m}{Z}^{-1}\delta ,$$

(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{x_i-x_0}{\Delta x}\right)}^{j-1},$$

(3)

in which x_i is the axial position of a stencil point, x₀ is the axial position of the point in which the derivative is demanded, and δ represents a vector such as:

$${\delta }_i=\{\begin{array}{cc}0& i\ne m+1\\ m!& i=m+1\end{array}.$$

(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 A_in is constructed as well. Both are used to evaluate the pure advection problem in the following form:

$$\frac{dT}{dt}=-u\left(AT+T_{in}{A}_{in}\right),$$

(5)

where T is the vector of localized fluid volume temperatures and T_in 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⁻⁶. 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).

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.

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 B_in 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}=-u\left(AT+T_{in}{A}_{in}\right)+D\left(BT+T_{in}{B}_{in}\right),$$

(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},$$

(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[\left(d_{1}B-A\right)T+\left(d_1{B}_{in}-A_{in}\right)T_{in}\right]=u\left(MT+M_{in}{T}_{in}\right),$$

(8)

where M and M_in 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 w_1,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 D_i 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).

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{array}{c}3.66\mathrm{if}Re\le 2300\\ 3.52\cdot U^4-45.148\cdot U^3+212.13\cdot U^2-427.45\cdot U^2+316.08\mathrm{if}2300\le Re\le 3100\\ \frac{f_r}{8}\cdot \frac{\left(Re-1000\right)\cdot Pr}{1+12.7\cdot {\left(\frac{f_r}{8}\right)}^{1/2}\cdot \left(P{r}^{2/3}-1\right)}\mathrm{if}Re<3000\end{array},$$

(9)

where Re is the Reynolds number, U = Re/1000, f_r 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},$$

(10)

where

$$\begin{array}{c}{\theta }_1={\left(-2.457\cdot \log \left({\left(\frac{7}{Re}\right)}^{0.9}+0.27\cdot \left(\frac{\epsilon }{d_1}\right)\right)\right)}^{16}\\ {\theta }_2={\left(\frac{37530}{Re}\right)}^{16}\end{array},$$

(11)

in which ε is the roughness of the inner wall and d₁ 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 \left(T_{sn,i,1}-T_i\right).$$

(12)

where T_sn,i,1 represents the inner wall (later first node of the solid region) temperatures of the solid shell, c_p 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:

$$P{e}_f=d_1\cdot u\cdot \frac{\rho \cdot c_p}{\lambda }\approx d_1\cdot u\cdot 6,29\cdot {10}^6\gg 1,$$

(13)

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

$$P{e}_s=d_1\cdot u\cdot \frac{{\rho }_s\cdot c_{p,s}}{{\lambda }_s}\approx d_1\cdot u\cdot 8\cdot {10}^4\gg 1.$$

(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{d{T}_{s,i}}{dt}\cdot {\rho }_s\cdot c_{ps}=\frac{1}{r}\frac{d}{dr}\left({\lambda }_s\cdot r\cdot \frac{d{T}_{s,i}}{dr}\right),$$

(15)

where r is the radial coordinate, T_s,i is the continuous temperature in the solid segment, ρ_s is the mass density, c_ps 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 \mathsf{\pi }\cdot \Delta x\cdot {\displaystyle \underset{r_1}{\overset{r_2}{\int }}r\cdot V\cdot \frac{d{T}_{s,i}}{dt}\cdot {\rho }_s\cdot c_{ps}\cdot dr}=2\cdot \mathsf{\pi }\cdot \Delta x\cdot \left[-{\displaystyle \underset{r_1}{\overset{r_2}{\int }}r\cdot \frac{dV}{dr}\cdot \frac{d{T}_{s,i}}{dr}\cdot {\lambda }_s\cdot dr}+{\left[r\cdot V\cdot {\lambda }_s\cdot \frac{d{T}_{s,i}}{dr}\right]}_{r_1}^{r_2}\right].$$

(16)

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

$$T_{se,i,j}\left(r\right)=\frac{\ln \left(r/R_j\right)\cdot \left(T_{sn,i,j+1}-T_{sn,i,j}\right)}{\ln \left(R_{j+1}/R_j\right)}+T_{sn,i,j}=T_{sn,i,j+1}\cdot \frac{\ln \left(r/R_j\right)}{\ln \left(R_{j+1}/R_j\right)}+T_{sn,i,j}\left(1-\frac{\ln \left(r/R_j\right)}{\ln \left(R_{j+1}/R_j\right)}\right),$$

(17)

where T_se,i,j are the temperature distributions in the cylindrical elements, T_sn,i,j are nodal temperatures, and R_j 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{d{T}_{se,i,j}}{dr}\right)=0\\ T_{se,i,j}\left(R_{j+1}\right)=T_{sn,i,j+1}\\ T_{se,i,j}\left(R_j\right)=T_{sn,i,j}\end{array},$$

(18)

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{array}{l}\{\begin{array}{c}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}\end{array},j=1\\ \{\begin{array}{c}\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}\end{array},2\le j\le n-1\\ \{\begin{array}{c}\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{array},j=n\end{array},$$

(19)

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:

$$C_{sn,j}=2\cdot \mathsf{\pi }\cdot \Delta x\cdot \left[{\rho }_{s,j-1}\cdot c_{ps,j-1}\cdot {\displaystyle \underset{R_{j-1}}{\overset{R_j}{\int }}r\cdot V_{sn,j}\cdot dr+{\rho }_{s,j}\cdot c_{ps,j}\cdot {\displaystyle \underset{R_j}{\overset{R_{j+1}}{\int }}r\cdot V_{sn,j}\cdot dr}}\right],$$

(20)

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

$$C_{sn,j}=\mathsf{\pi }\cdot \Delta x\cdot \left[{\rho }_{s,j-1}\cdot c_{ps,j-1}\cdot \left(\frac{R_{j-1}^2-R_j^2}{2\cdot \ln \left(R_j/R_{j-1}\right)}+R_j^2\right)-{\rho }_{s,j}\cdot c_{ps,j}\cdot \left(\frac{R_j^2-R_{j+1}^2}{2\cdot \ln \left(R_{j+1}/R_j\right)}+R_j^2\right)\right],$$

(21)

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_{sn,i,j}=2\cdot \mathsf{\pi }\cdot \Delta x\cdot \left[T_{sn,i,j-1}\cdot h_{sn,j-1,j}+T_{sn,i,j}\cdot h_{sn,j,j}+T_{sn,i,j+1}\cdot h_{sn,j+1,j}\right],$$

(22)

where

$$\begin{array}{c}{h}_{sn,j-1,j}=-{\displaystyle \underset{R_{j-1}}{\overset{R_j}{\int }}r\cdot {\lambda }_{s,j-1}\cdot \frac{d{V}_{sn,j-1}}{dr}\frac{d{V}_{sn,j}}{dr}dr}=\frac{{\lambda }_{s,j-1}}{\ln \left(R_j/R_{j-1}\right)}\\ h_{sn,j,j}=-{\displaystyle \underset{R_{j-1}}{\overset{R_j}{\int }}r\cdot {\lambda }_{s,j-1}\cdot {\left(\frac{d{V}_{sn,j}}{dr}\right)}^{2}dr-{\displaystyle \underset{R_j}{\overset{R_{j+1}}{\int }}r\cdot {\lambda }_{s,j}\cdot {\left(\frac{d{V}_{sn,j}}{dr}\right)}^{2}dr}}=\frac{-{\lambda }_{s,j-1}}{\ln \left(R_j/R_{j-1}\right)}+\frac{-{\lambda }_{s,j}}{\ln \left(R_{j+1}/R_j\right)}.\\ h_{sn,j,j}=-{\displaystyle \underset{R_j}{\overset{R_{j+1}}{\int }}r\cdot {\lambda }_{s,j}\cdot \frac{d{V}_{sn,j+1}}{dr}\frac{d{V}_{sn,j}}{dr}dr}=\frac{{\lambda }_{s,j}}{\ln \left(R_{j+1}/R_j\right)}\end{array}$$

(23)

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_{sn,i,1}=2\cdot \mathsf{\pi }\cdot \Delta x\cdot \left[T_{sn,i,1}\cdot h_{sn,1,1}+T_{sn,i,2}\cdot h_{sn,2,1}\right]+\pi \cdot \Delta x\cdot Nu\cdot \lambda \cdot \left(T_i-T_{sn,i,1}\right),$$

(24)

where

$$\begin{array}{c}{h}_{sn,1,1}=-{\displaystyle \underset{R_1}{\overset{R_2}{\int }}r\cdot {\lambda }_{s,1}\cdot {\left(\frac{d{V}_{sn,1}}{dr}\right)}^{2}dr}=\frac{-{\lambda }_{s,1}}{\ln \left(R_2/R_1\right)}\\ h_{sn,2,1}=-{\displaystyle \underset{R_1}{\overset{R_2}{\int }}r\cdot {\lambda }_{s,1}\cdot \frac{d{V}_{sn,2}}{dr}\frac{d{V}_{sn,1}}{dr}dr}=\frac{{\lambda }_{s,1}}{\ln \left(R_2/R_1\right)}\end{array},$$

(25)

and similarly for the outer boundary node:

$$H_{sn,i,n}=2\cdot \mathsf{\pi }\cdot \Delta x\cdot \left[T_{sn,i,n-1}\cdot h_{sn,n-1,n}+T_{sn,i,n}\cdot h_{sn,n,n}\right]+{\alpha }_{out}\cdot 2\cdot \pi \cdot R_n\cdot \Delta x\cdot \left(T_{air}-T_{sn,i,n}\right),$$

(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, T_air is the ambient temperature and:

$$\begin{array}{c}{h}_{sn,n-1,n}=-{\displaystyle \underset{R_{n-1}}{\overset{R_n}{\int }}r\cdot {\lambda }_{s,n-1}\cdot \frac{d{V}_{sn,n-1}}{dr}\frac{d{V}_{sn,n}}{dr}dr}=\frac{{\lambda }_{s,n-1}}{\ln \left(R_n/R_{n-1}\right)}\\ h_{sn,n,n}=-{\displaystyle \underset{R_{n-1}}{\overset{R_n}{\int }}r\cdot {\lambda }_{s,n-1}\cdot {\left(\frac{d{V}_{sn,n}}{dr}\right)}^{2}dr}=\frac{-{\lambda }_{s,n-1}}{\ln \left(R_n/R_{n-1}\right)}\end{array}.$$

(27)

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

$$\frac{d{T}_{sn,i,j}}{dt}=\frac{H_{sn,i,j}}{C_{sn,j}},$$

(28)

Values in matrix h_sn,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 \mathsf{\pi }\cdot \Delta x\cdot \frac{h_{sn,i,j}}{C_{sn,j}},$$

(29)

$$k_{in,i,j}=\{\begin{array}{c}j=1,\mathsf{\pi }\cdot \Delta x\cdot Nu\cdot \lambda \cdot \left(T_i-T_{sn,i,1}\right)/C_{sn,1}\\ 2\le j\le n-1,0\\ j=n,{\alpha }_{out}\cdot 2\cdot \mathsf{\pi }\cdot R_n\cdot \Delta x\cdot \left(T_{air}-T_{sn,i,n}\right)/C_{sn,n}\end{array}.$$

(30)

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

$$\frac{d{T}_{sn,i}}{dt}=k\cdot T_{sn,i}+k_{in,i}.$$

(31)

where T_sn,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 k_in,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{d{T}_{sn}}{dt}=K\cdot T_{sn}+K_{in}.$$

(32)

where T_sn is the column vector filled with column vectors T_sn,i, K is the banded square matrix filled with k on its diagonal, and K_in is the time-dependent boundary vector filed with column vectors k_sn,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 property material of the steel circular wall (inner pipe diameter is 0.6 m).

Property	Case 1 Accumulation	Case 2 No Accumulation
Layer thickness (mm)	100	4
Mass density (kg/m3)	8030	8030
Heat conductivity (W/mK)	16.27	16.27
Specific heat capacity (J/kgK)	502.45	502.45

). 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.

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 Figure 8, Figure 9, Figure 10 and Figure 11 with Figure 12, Figure 13, Figure 14 and Figure 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 Figure 8, Figure 9, Figure 10 and Figure 11). The whole event of heat recharging is used (7200 s). A tolerance of 10⁻¹² is chosen to ensure that the time stepping method can always take advantage of a superior spatial discretization.

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 (Figure 13 and Figure 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 10⁻⁶, 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 Figure 16 and Figure 17). The computational time can be mapped directly (see Figure 18).

Similar tests for the second case clearly demonstrate the independence on the radial number of subdivisions (see Figure 19 and Figure 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.

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 Figure 23 and Figure 24) and STAR-CCM+ (fully 3D in Figure 25 and Figure 26). The tolerance of the new model was lowered to 10⁻⁵ 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.

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. Figure 27 and Figure 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.

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 Figure 23, Figure 24, Figure 25 and Figure 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

• The shape function used in the radial model as well as the use of analytic description of heat exchange between the fluid and the inner wall surface ensures that the solution of the new model agree with the analytic solution completely even for the roughest radial discretization (also visible through fast stabilization in Figure 22 after the sudden temperature change).
• The new model allows simple modeling of circular pipes with an arbitrary number of shells and their material properties (e.g., a steel wall followed by insulation and a plastic shell). All shells then consider radial dynamic effects based on the newly derived radial model, and OpenModelica’s marching algorithm automatically adapts the time step to ensure its numerical stability.
• Although elements growing in size in the radial direction better capture the sudden changes of temperature, the necessity of smaller time steps to ensure the stability of the marching algorithm completely negates this advantage (compare the CPU time in Figure 7) and makes the uniform mesh more preferable.
• With the preferred use of the uniform mesh in the radial model, the thick wall has to be modeled with higher number of elements, while the thin wall is captured well with very few or even single element (compare error estimations for grow factor 1 in Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13, Figure 14 and Figure 15 or the error maps in Figure 16 and Figure 19).
• The use of the method of lines and appropriate algorithm with adaptable time step makes the precision easily tradable for computation time, without the necessity to focus on the numerical stability conditions of the composite model (trading the capturing of the transients for CPU time, the steady-states will be captured accurately even with rough radial meshes).
• There are regions in which a sole improvement of the mesh in one direction will not result in superior solution (e.g., see the region of 100 axial and 20 radial elements in Figure 16, where sole improvement of radial discretization does not lower errors). There is, however, an ideal direction for the discretization in each case.
• Figure 22 shows how does strong thermal inertia of the wall affect the outlet temperature in idealized case (the effect lasts approximately for half an hour for the thick steel wall).
• Maximal deviation in outlet temperature in comparison with Fluent (2D) is about 3 K.
• Maximal deviation in outlet temperature in comparison with STAR-CCM+ (3D) is about 6 K for the case without gravity and 4 K for the case that includes gravity.
• During all simulations, positive weights in the definition of axial diffusion are deployed, which shows that turbulent axial mixing plays a role in thermal wave propagation.
• Presence of gravity in fully horizontal pipes acts as additional axial diffusion from the viewpoint of cross-sectional temperature average (comparison of the two curves from STAR-CCM+ in Figure 26).
• Relative simulation time of the new model for a single 100 m long pipe with strong accumulation is approximately in order of 3.75 × 10⁻⁵.
• A simple quasi-static evaluation of the pressure losses (the friction coefficient is already being evaluated during thermal simulation) would allow the use of the new model for simulation of whole grids with automatic mass flows evaluation in each branch (e.g., OpenModelica connector definition includes Kirchhoff laws by default).

5. Conclusions

This paper deals with transients in circular pipes, where the accumulation of heat into surrounding material has been considered in detail. Model is a combination of the analytic description of a heat transfer between the fluid and the solid wall, simple advection scheme, turbulent diffusion, and a finite element scheme dealing with the circular wall in which the shape functions are based on the steady-state solution. Method of lines is used for time integration. The new model represents an effective numerical solution to the thermal wave propagation problem, which includes advanced evaluation of thermal inertia of the pipe’s wall and the axial diffusion that is caused by turbulence. The diffusion coefficient evaluation is mainly proportional to flow velocity.

Overall, it can be said that strong accumulation affects the thermal wave propagation for a stretched period of time in comparison with the non-accumulation case. The effect of the re-accumulation is not only dependent on the pure amount of re-accumulated heat, but also on the particular way this event happens. In other words, to capture the event accurately is to focus on the proper evaluation of the innermost surface temperature of the pipe’s wall because it directly affects the rate in which the temperature of the fluid changes. The discretization of both axial and radial direction must go hand in hand to improve the solution. During comparison with the fine models from Fluent and STAR-CCM+, deployment of positive weights in the diffusion coefficient evaluation are always necessary to match the simulated data. This means that turbulent diffusion plays a role in temperature wave propagation. Presence of physical diffusion also reduces unphysical oscillation on sharp interfaces. Gravity is shown to resemble the behavior of an additional axial diffusion as well. The new model resembles the transient behavior of the accumulating pipe wall with highly comparable results to the equivalent, very fine models from Fluent and STAR-CCM+, which should represent the ideal physics-based modeling. The validation with a branched real system needs to be done in the future.

The composite model had run faster in OpenModelica than in Julia possibly because of a non-optimized jacobian evaluation. It might be possible to speed up the Julia execution using function definition macros from Differential Equations package, which could be able to derive symbolic jacobian automatically. The model can be extended by pressure loss evaluation, which would allow simulation of whole DH grids based on Kirchhoff laws. Regularization of mass flows would be necessary here to avoid singularities (division by zero) in the nodes.

The future work should aim at a self-identification algorithm using inlet/outlet temperature signals from real installations to calibrate the appropriate weights in the model. It may be possible to generalize the approach for more complex configurations, such as twin-pipes, buried pipes, and so on. An optimization algorithm can be used to find the matrix coefficient describing the radial heat dynamics. A background from machine learning may be utilized for such an endeavor. The scaling of the new model, when used for simulation of whole grids with automatic mass flow evaluation in each branch based on the pressure drops between individual nodes, is also part of the future work.