New Dimensionless Correlations for the Evaluation of the Thermal Resistances of a District Heating Twin Pipe System

: The heat losses from pre-insulated double-pipe district heating (DH) systems buried in a homogeneous soil are studied numerically. The study is conducted using the diameter of the pipes and their distance, the size of the insulation, the thermal conductivity ratio between the insulation and the soil, as well as the burial depth of the double-pipe system, as controlling parameters. A computational code based on a control-volume formulation of the ﬁnite-difference method has been developed using the open-source framework OpenFOAM with the purpose to compute the heat transfer rate across adjacent solid regions with different thermophysical properties. The main scopes of the study are: (a) to investigate in what measure the geometry and the relative position of the warm and cold pipes, as well as the temperature imbalance, the burial depth and the physical properties of both the insulation and the soil, affect the heat losses; (b) to analyze the existence of an optimal conﬁguration of the DH system by the thermal resistance enhancement viewpoint; and (c) to develop accurate correlating equations for the evaluation of the thermal resistance existing between each pipe and its surroundings, useful for practical thermal engineering applications.


Introduction
The development of a network of pipes connecting buildings for heat supply by renewable source and/or heat recovery from industrial surplus has rapidly grown in Europe in the last ten years, as reflected by the thousands of district heating (DH) systems covering about 9% of the European thermal energy demand [1]. On the other hand, with the increasing introduction of such systems, the national governments need specific requirements for reducing the environmental impact of the future smart energy grids, and to be more competitive in developing energy saving strategies. In this scenario, one of the main challenges is the transformation of the current DH systems into the fourth generation district heating (4GDH), with the implementation of renewable energy sources with lower and more flexible distribution temperatures, and the replacement of pressurized hot water and steam with low-enthalpy fluids [2].
The advantages in using lower temperatures for heat distribution over long pipes are primarily in the reduction of the thermal losses between the pipes and the ground, which increases the DH thermal efficiency. On the other hand, even with supply temperatures near 50 • C, the typical energy losses can reach 20% of the total heat generation, mainly depending on the weather conditions, the insulation thermal conductivity and thickness, as well as the depth and orientation of the pipes with respect to the ground surface [3]. Indeed, in the last few years, many efforts have been made to increase the thermal resistance around the pipes, which has led to the pre-insulated double-pipe arrangements currently in use.
In pre-insulated double-pipe systems, the total thermal losses are primarily affected by two opposite effects arising from the relative position of the warm and cold pipes within the insulation: the first effect, which tends to increase the heat transfer toward the external environment, prevails when the distance between the pipes increases, thus reducing the thermal resistance between each pipe and the insulation casing; the second effect, which tends to reduce the total thermal losses, prevails when the pipes are at a small distance and the heat transfer rate from the supply pipe to the return pipe increases. The situation is such that, to minimize the heat losses from the supply pipe, an optimum distance between the pipes does exist.
The correlations typically used for the calculation of the effective thermal resistance of one or two buried pipes, based on the analytical solutions developed by Carslaw and Jaeger [4], are those developed by Wallentén [5], Claesson and Bennet [6], and Claesson and Javed [7]. However, the proposed equations, derived by the complex mathematical formulation of the multipole method, provide accurate solutions only in the hypothesis of a combination of cylindrical sources located inside a cylindrical domain with a uniform temperature condition imposed at its boundary, which is not strictly verified in real situations.
The existence of an optimal distance between two pipes at different temperatures embedded in a cylindrical insulation has been investigated experimentally by Wonsyld et al. [8], who found that the heat losses from the supply pipe decrease when the supply pipe is located at the center of the pipeline, whereas the return pipe is located closer to the external surface of the insulation. Later, Bøhm and Kristjansson [9] calculated both analytically and numerically the heat losses occurring in three sample-specific cases represented by a pair of uninsulated pipes, a pair of pipes embedded in a cylindrical insulation and a pair of pipes embedded in an egg-shaped insulation, with the only aim being to demonstrate the better performance of the egg-shaped configuration. Dalla Rosa et al. [10] carried out FEM simulations for buried pre-insulated double pipes in the hypothesis that the insulation thermal conductivity was temperature dependent, again showing that the thermal resistance between the supply pipe and the ground surface is larger for asymmetrical configurations, that is, when the pipes are located at different distances from the center of the insulation, although such a result is affected by the imposition of a radial thermal field at a close distance from the insulation casing. More recently, Danielewicz et al. [11] performed a three-dimensional simulation of the temperature field around two buried pipes, inside which the forced convection of the flowing water was taken into account by means of the Reynolds-averaged Navier-Stokes (RANS) equations. The predicted temperature distributions along the pipes due to the mutual interactions occurring between the hot fluid, the cold fluid and the ground were compared with some available experimental data, showing a good agreement. However, although the numerical model was undoubtedly able to capture the local heat losses through the soil, the applicability of the results obtained is limited to the specific case analyzed. The temperature distributions and the associated thermal losses along double-pipe systems have also been determined by van der Heijde et al. [12], who applied the aforementioned correlations of Wallentén in conjunction with the energy balance equations for counter flow heat exchangers, despite that such an approximate approach gives rise to non-negligible discrepancies when the mutual influence between the supply and return pipes are considered. Ocłøn et al. [13] used a finite-element method in a rectangular integration domain to simulate the temperature distribution in a district heating pipeline and the surrounding ground, for various types of insulation, both in winter and summertime climatic conditions. Finally, in the study conducted by Krawczyk and Teleszewski [14], the heat losses from a double-pipe system embedded in either a circular or an elliptic insulation were calculated numerically. Besides the conclusion that the elliptic geometry reduces the overall heat losses, the main result is the existence of an optimal distance between the supply and return pipes allowing for the smallest heat loss from the supply pipe and, at the same time, a limited heat transfer rate from the warm supply pipe to the cold return pipe. Other studies with a bearing on the subject are those performed by Schneider [15], Bøhm [16], Afshan and Pettinger [17] and Teleszewski et al. [18]. In conclusion, it seems worth pointing out that, leaving aside the approximate equations derived by the multipole method [5][6][7], no other correlation was developed, based either on numerical data or experimental results, to more accurately predict the thermal resistance existing between each tube of the insulated pipeline and its surroundings.
Framed in this general background, a comprehensive numerical study on pre-insulated double-pipe DH systems buried in a homogeneous soil is executed by solving the steadystate equation of heat conduction by way of the finite-difference method. The study is conducted using the diameter and location of the pipes, the size of the insulation, the temperature ratio between the warm and cold pipes, the thermal conductivity ratio between the insulation and the soil, as well as the burial depth of the double-pipe system, as controlling parameters. The main scopes of the study are: (a) to investigate in what measure the geometry and the relative position of the warm and cold pipes, as well as the temperature imbalance, the burial depth and the physical properties of both the insulation and the soil, affect the heat losses; (b) to analyze the existence of an optimal configuration of the DH system by the thermal resistance enhancement viewpoint; and (c) to develop accurate correlating equations for the evaluation of the thermal resistance existing between each pipe and its surroundings, useful for practical engineering applications.

Mathematical Formulation
A pre-insulated double-pipe DH system is considered. The pipeline consists of a pair of horizontal circular tubes of a diameter of d, embedded in a cylindrical insulation having a diameter of w, buried at a depth of h beneath the ground surface, which represents the upper boundary of the rectangular integration domain, as sketched in Figure 1, where the reference polar and Cartesian coordinate systems for the pipes, the insulation and the soil domain are also displayed. The pipes, set side by side at a center-to-center distance of x p , are kept at the uniform temperatures of t s and t r , in which the subscripts "s" and "r" denote the supply pipe and the return pipe, respectively. As far as the boundaries of the integration domain are concerned, the ground surface is kept at a uniform temperature of t g lower than both t s and t r , while the other pseudo-boundaries are assumed to be adiabatic. The resulting thermal field is considered to be steady and two-dimensional, whereas the thermal conductivities of the insulation and the ground, k i and k g , respectively, are assumed to be constant. Under these hypotheses, the steady-state governing equation of conduction heat transfer expressed in dimensionless form reduces to where T is the dimensionless temperature excess over the uniform temperature of the ground surface normalized by the temperature difference (t s − t g ) The related boundary conditions are: (a) T s = 1 and T g = 0 at the outer surface of the supply pipe and at the ground surface, respectively, (b) 0 < T r ≤ 1 at the outer surface of the return pipe, (c) ∂T/∂Y = 0 at the bottom horizontal pseudo-boundary line and (d) ∂T/∂X = 0 at both of the vertical pseudo-boundary lines, where X and Y are the dimensionless horizontal and vertical Cartesian coordinates normalized by the tube diameter d, respectively. Additionally, the heat flux and temperature continuity at the interface between the insulation and the surrounding ground are also imposed, which means k i · (∂T/∂n) i = k g · (∂T/∂n) g and T i = T g where n denotes the normal to the interface, while the subscripts "i" and "g" denote the insulation and the ground, respectively.
Other dimensionless parameters which enter into this study are: (a) The dimensionless diameter of the insulation The dimensionless center-to-center pipe spacing (c) The dimensionless burial depth of the double-pipe system (d) The thermal conductivity ratio between the insulation and the ground Other interesting configurations, such as the vertical alignment of the pipes or asymmetrical locations of the supply and return pipes with respect to the center of the insulation, will be the object of future investigations.

Computational Procedure
The governing equation, in conjunction with the boundary conditions stated earlier, is solved by means of a control-volume formulation of the finite-difference method using the open-source framework OpenFOAM [19]. A new solver, derived from the combination of the basic solvers named scalarTransportFoam and chtMultiRegionFoam, has been developed with the purpose to compute the heat transfer across adjacent solid regions with different thermophysical properties. According to the geometry of the system, fine cylindrical polar grids are used in close proximity of both the outer surface of each tube and the inner and outer surface of the insulation, whereas the remainder of the interior of the insulation and the surrounding ground at a short distance are filled with linear triangular elements. Conversely, a Cartesian grid is used to discretize the integration domain at a larger distance from the pipeline, as depicted in Figure 2. Starting from a uniform temperature distribution across the integration domain, the discretized governing equation is solved iteratively by way of the conjugate gradient algorithm using a second-order central difference scheme. A standard under-relaxation technique is enforced in all of the steps of the computational procedure to ensure an adequate convergence. The solution is considered to be converged when the temperature change at any grid node between two consecutive iterations is smaller than 10 −6 , and the difference between the incoming and outgoing heat transfer rates, at the pipes and the ground surface, is smaller than 10 −2 .

Pipe
Once a converged solution is achieved, the dimensionless heat transfer rates q s and q r at the outer surfaces of the supply and return pipes, normalized by the product of the temperature difference (t s − t g ) and the insulation thermal conductivity k i , are calculated as: On the other hand, the dimensionless heat transfer rates q s and q r can also be expressed as: where ρ pg and ρ pp are the dimensionless thermal resistances existing between each pipe and the ground surface and between the two pipes, respectively, normalized by the reciprocal of the insulation thermal conductivity k i . Therefore, once the temperature of the return pipe is assumed to be the same as that of the supply pipe, i.e., T r = T s = 1, the combination of Equations (7) and (9) allows for the calculation of the value of ρ pg . Clearly, the same result can be achieved by using Equations (8) and (10). Subsequently, once the value of ρ pg is replaced into Equations (9) and (10) and the temperature of the return pipe is assumed to be lower that of the supply pipe, e.g., T r = 0.5, the combination of Equations (7) and (9) allows for the calculation of the value of ρ pp . Of course, the same result can be reached whatever value lower than unity is assumed for T r , as well as using Equations (8) and (10) instead of Equations (7) and (9).
Numerical tests on the dependence of the results on the mesh spacing, and on the extent of the whole computational domain, have been methodically performed for several combinations of the main controlling parameters. In particular, the optimal number of mesh elements and the optimal positions of the pseudo-boundary lines used for computations are such that further grid refinements or domain expansions do not produce changes larger than 1% in both the heat transfer rates and the predicted temperature field.
The selected results of the grid sensitivity analysis obtained for different values of both K and T r are listed in Table 1, in which the effects of the X-wise and Y-wise extents of the computational domain, as well as the number of mesh elements, on the value of the dimensionless thermal resistances ρ pg and ρ pp have been calculated for successive enlargements of the integration domain and mesh refinements. The typical numbers of the nodal points used for simulations lie in the ranges between 80,000 and 200,000 elements, whilst the X-wise and Y-wise extents of the whole integration domain are up to 20 times and 10 times the burial depth H, respectively.
Finally, with the scope to validate the numerical code used for the present study, two different tests have been executed. In the first test, the solutions obtained for the dimensionless thermal resistance ρ pg existing between each pipe of a pair of buried pipes kept at the same temperature and the ground surface have been compared with the corresponding values derived using the analytical solutions provided by Carlslaw and Jaeger [4]. In the second test, the solutions obtained for the dimensional heat transfer rates for unit length occurring between each pipe of a pre-insulated double-pipe system and the insulation casing have been compared with the numerical data of Krawczyk and Teleszewski [14]. In both cases, a good degree of agreement between our numerical solutions and the literature data was achieved, as shown in Tables 2 and 3.

Results
The numerical simulations are performed for different values of (a) the dimensionless diameter of the insulation, W, in the range between 2 and 6, (b) the dimensionless centerto-center pipe spacing, L, in the range between 1.2 and 3.6, (c) the dimensionless burial depth, H, in the range between 1 and 10, (d) the ratio between the thermal conductivities of the insulating material and the ground, K, in the range between 0.01 and 0.5 and (e) the dimensionless temperature of the return pipe, T r , in the range between 0 and 1.
The selected local results are displayed in Figures 3 and 4, in which equispaced isotherm contours are plotted for H = 6, W = 4 and L = 1.3, 1.75 and 2.5. As expected, the temperature field in the ground is featured by an eccentric distribution with a higher isotherm concentration in the region between the ground surface and the pre-insulated double pipe. On the other hand, the highest temperature gradients are localized within the insulation, in which the distance between the pipes plays the most important role in determining the heat transfer rate from the supply pipe to the return pipe, and from each pipe to the ground surface. In fact, for small center-to-center spacings, most of the heat loss from the supply pipe is "gained" by the return pipe; conversely, when the pipe spacing is large enough, the higher insulation thickness between the pipes results in a growth of the amount of heat transferred to the ground surface. Therefore, the thermal resistance existing between each pipe and the ground surface is affected by the heat transfer rate occurring from the supply pipe to the return pipe, which means that ρ pg increases as L is decreased. This is clearly illustrated in Figures 5 and 6, in which the heat flux at the supply and return pipe surfaces are plotted for L = 1.3, 1.75 and 2.5. Additionally, ρ pg is expected to increase as the burial depth H is increased, as reported in Figure 7. In parallel, the thermal resistance ρ pp existing between the supply and return pipes increases as L is increased and H is decreased, as displayed in Figure 8. Hence, an optimal dimensionless pipe spacing of L opt for the minimum heat loss from the supply pipe for a prescribed temperature imbalance assigned between the pipes has to be expected, as displayed in Figure 9, where a number of distributions of the dimensionless heat transfer rate from the supply pipe q s are plotted versus L for different values of T r , showing that the optimal pipe spacing decreases as the temperature of the return pipe increases.     The effects of W, H and K on the heat loss from the supply pipe are then pointed out in Figures 10-12 in which, for any independent variable, a number of distributions of q s are plotted versus L using the variable itself as a parameter. It is apparent that the minimum heat loss increases as the insulation diameter decreases, whereas all of the other controlling parameters have more limited effects, at least within their investigated ranges of variability. Furthermore, it can be noticed that the optimal pipe spacing increases as the insulation diameter and the burial depth are increased.
Based on the results obtained, the dimensionless thermal resistances ρ pg and ρ pp can be expressed as a function of all of the considered dimensionless controlling parameters through the following empirical correlations derived by a multiple regression method: ρ pp = 0.016 · L 4.26 H −0.155 W −2.33 K −0.66 + 0.72 · L 0.867 − 0.5 · H 0.086 (12) The standard deviation of the errors of Equations (11) and (12) are 3.2% and 3.5%, respectively, with a percent range of error of ±5% and a 95% level of confidence, as shown in Figures 13 and 14.   Notice that the results obtained by the employment of these correlations differ by less than 5% from the numerical solutions computed by Bøhm and Kristjansson [9] for a buried pre-insulated double-pipe system, as reported in Table 4, which demonstrates that Equations (11) and (12) can safely be used to predict the heat losses from underground twin pipes in a wide range of application for district heating systems. Table 4. Comparison of the present solutions with the numerical data of Bøhm and Kristjansson [9] for a buried pre-insulated double-pipe system.  Finally, the dimensionless optimal center-to-center pipe spacing L opt can be correlated to the controlling parameters by the following empirical equation: with a 3.1% standard deviation of errors, with a percent range of error of ±5% and a 98% level of confidence, as shown in Figure 15.

Conclusions
A comprehensive numerical study on pre-insulated double-pipe systems buried in a homogeneous soil has been executed by solving the steady-state equation of heat conduction by way of the finite-difference method. In particular, the dimensionless thermal resistances between each pipe and the ground surface, and between the two pipes, have been calculated with the scope to obtain the global heat losses occurring through the soil. Simulations have been performed using the center-to-center distance between the pipes and the dimensionless burial depth of the pre-insulated double-pipe system, both normalized by the pipe diameter, as well as the ratio between the thermal conductivities of the soil and the insulation, as independent variables.
The main results obtained in the present study may be summarized as follows: (a) The existence of an optimal distance between the pipes which minimizes the heat losses from the supply pipe has been found for any investigated set of the controlling parameters; (b) The optimal pipe spacing increases with increasing both the insulation diameter W and the burial depth H, whereas the optimal pipe spacing decreases as the temperature T r of the return pipe increases; (c) The minimum heat loss increases as the insulation diameter decreases; (d) A set of correlations to predict the heat losses from underground double-pipe systems, as well as the optimal pipe spacing, has been developed for wide ranges of the main controlling parameters.