Relation Between Mean Fluid Temperature and Outlet Temperature for Single U-Tube Boreholes

Ground-coupled heat pump (GCHP) systems usually utilize buried vertical heat exchangers, named borehole heat exchangers (BHEs). The accurate design or simulation of a GCHP system requires the calculation of the time-dependent outlet temperature from the BHEs, Tout. However, the most widely employed BHE simulation models yield the time evolution either of the mean temperature of the BHE-ground surface, Tsm, or of that of the fluid, Tfm. In transient regime, it is not easy to relate Tout to either Tsm or Tfm. In this paper we determine, through 3D finite element simulations, simple expressions of a dimensionless coefficient φ allowing the calculation of Tout by means of a simulation model that yields Tfm. These expressions hold for single U-tube BHEs, both in quasi-steady and in unsteady working conditions. We validate our 3D simulation code by comparison with an analytical BHE model. Then, we present applications of our expressions of φ to calculate the time-dependent values of Tout through a BHE model that yields those of Tfm. Finally, we show that the values of φ in quasi-steady working conditions can be used for a simple calculation of the effective borehole thermal resistance.


Introduction
The diffusion of ground-source heat pumps, and in particular, of ground-coupled heat pumps (GCHPs), is rapidly increasing. Indeed, GCHPs are a very efficient technology for the climatization of buildings and the production of domestic hot water [1]. The performance of these systems has been analyzed by experiments and by simulation tools [2][3][4][5][6][7]. GCHPs usually utilize buried vertical heat exchangers, named borehole heat exchangers (BHEs), that are mostly composed of a single or double polyethylene U-tube, installed in a drilled hole subsequently filled with a grouting material. A BHE has a length that commonly ranges from 50 to 200 m, and a diameter of about 15 cm. Improving the BHE-field design is an important way to enhance the efficiency of GCHP systems.
Several models for a BHE-field design have been proposed [8][9][10][11][12]. Most methods adopt dimensionless factors of thermal-response, also named g-functions. The g-function of a borefield yields the time-dependent dimensionless mean temperature of the boundary surface of the borefield caused by a constant total power released by the BHEs. Different methods to obtain the g-function of a borefield have been developed, either assuming that the heat flux released by the BHEs is uniform along the surface between BHEs and soil [13][14][15][16][17], or assuming that the total heat flux released by the BHEs is constant in time and the boundary between the ground and the borefield has a uniform temperature [18][19][20][21]. More accurate thermal response factors have also been determined, by assuming time-constant total heat flux and same inlet temperature in all the BHEs, and considering the influence of the borehole thermal resistance on the borefield surface temperature [22,23].
The g-functions or thermal response factors determined by the methods cited above are not precise in the short term, i.e., during the first one or two hours of operation, because the models employed do not consider accurately the thermal inertia of the BHE. Moreover, they yield the surface averaged time-dependent temperature of a borefield, T sm , but not the time-dependent mean temperature of the fluid, T fm . The latter is then determined by the following equation: where R b is the BHE thermal resistance per unit length and q lm is the mean linear thermal power exchanged between the fluid and the surrounding ground. Equation (1) is correct when the heat transfer in the borehole is quasi-steady, but is not precise in the short term and yields an unphysical jump of T fm when q lm changes from zero to a given constant value.
In order to predict accurately the time evolution of T fm during the first hours of operation, several researchers developed simulation models that are accurate even in the short term.
De Carli et al. [24] and Zarrella et al. [25] developed a Capacity Resistance Model (CaRM) of BHE suitable for the short-term analysis and employable also in the long term. Quaggiotto el al. [26] applied the CaRM model proposed in [25] for a numerical and experimental comparison between coaxial and double U-tube BHEs. Other resistance and capacity models were presented by Bauer et al. [27] and by Pasquier and Marcotte [28]. Ruiz-Calvo et al. [29] proposed to separate the short-term and the long-term simulation and developed a short-term model based on that by Bauer et al. [27]. Li and Lai [30,31] developed a 2D analytical BHE model where the tubes are schematized as infinite line sources that supply a uniform and constant linear heat flux. Zhang et al. [32] presented a transient quasi-3D line source model that introduces the concept of transient borehole thermal resistance and gives a full-time-scale thermal response.
Beier and Smith [33] proposed an analytical cylindrical BHE model, where the borehole is represented by a grout annulus with external radius equal to that of the borehole, and internal radius such that the grout annulus has a thermal resistance equal to that of the borehole.
Xu and Spitler [34] developed a numerical model that approximates the real BHE structure with several concentric cylinders, that include a fluid layer, an equivalent convective-resistance layer, a tube layer, and a grout layer. Man et al. [35] presented the analytical solution for a simple cylindrical BHE model, where the thermal properties of the BHE materials coincide with those of the ground and the thermal power is supplied by a generating cylindrical surface that represents the fluid. The solution is given both for the 1D scheme, that neglects the heat conduction along the BHE axis, and for the 2D axisymmetric scheme, that takes into account the finite length of the BHE.
Bandyopadhyay et al. [36] found an analytical solution for a borehole model composed of a high-conductivity solid cylinder with heat generation, representing the fluid, and a grout layer. The solution was given in the Laplace transformed domain, and the authors employed a numerical inversion to determine the thermal response in the time domain. Javed and Claesson [37] developed a complete analytical solution of the borehole model employed by Bandyopadhyay et al. [36]. Claesson and Javed [38] proposed an analytical method to determine the thermal response of a borehole, valid both in the short term and in the long term, by coupling the short-term model presented in Javed and Claesson [37] to a long-term model based on the finite line-source solution.
Beier [39] developed an analytical BHE model that approximates the U-tube as two half tubes, and considers both the fluid and the grout thermal capacity. The analytical solution is given in the Laplace transformed domain and is inverted by using the Stehfest algorithm. Lamarche [40] developed an analytical cylindrical borehole model in which the borehole is composed of a solid cylinder subjected to heat generation, representing the fluid, surrounded by cylindrical layers representing the polyethylene pipes, and the grout. Naldi and Zanchini [41] proposed a numerical BHE model (OMEC) composed of a homogeneous equivalent cylinder having an internal heat-generating surface and thermal properties suitable to reproduce both the thermal resistance and the heat capacity of the BHE. By means of that Energies 2020, 13, 828 3 of 23 model, the authors determined the full-time-scale evolution of T fm for a borefield with BHEs having equal inlet temperatures.
Most of the short-term or full-time-scale simulation models cited above yield directly the time evolution of T fm , for a borefield subjected to a time constant heat load, without employing Equation (1). However, they do not yield directly the time evolution of the outlet fluid temperature, T out , that is needed for the simulation of the heat pump. Therefore, it is interesting to complete these models by providing relations between T fm and T out .
Beier and Spitler [42] presented a method that allows determining the dimensionless factor f, given by: where T in is the inlet fluid temperature. Equation (2) and the energy balance equation: where .
Q is the thermal power supplied by the heat pump to the borehole fluid, . m and c p are the mass flow rate and specific heat capacity of the fluid, yield: Through 3D finite element simulations and best fit of simulation results, Zanchini and Jahanbin [43] determined simple correlations to evaluate the dimensionless coefficient ϕ defined as: where .
V is the volume flow rate, . V 0 is a reference value of . V, namely 12 L per minute (L/min), and T ave = (T in + T out )/2. The correlations reported in [43] apply to double U-tube boreholes. Equation (5) and the balance Equation (3) yield: In the present paper, new correlations to determine ϕ are provided, for single U tube BHEs, through the best fit of the results of 3D numerical simulations. These correlations can be employed to obtain an accurate evaluation of the time evolution of T out by means of a BHE simulation code that yields the time evolution of T fm . An example of this use is reported, in the case of constant flow rate and constant power supplied to the ground, by evaluating the time evolution of T fm through the simple analytical BHE model proposed by Man et al. [35], and that of T out through Equation (6). Then it is shown that our correlations for ϕ can be applied to determine accurately the time evolution of T out from that of T fm even in the simulation of BHE fields subjected to a time dependent heat load. Finally, it is shown that the correlation for ϕ valid for the quasi-stationary regime can be employed for an immediate calculation of the effective BHE thermal resistance.

Simulation Cases and Method
We considered nine BHE geometries, with BHE diameter D b = 152 mm, shank spacing d = 84, 94 and 104 mm, length L = 50, 100 and 200 m, pipes with internal diameter D pi = 32.6 mm and external diameter D pe = 40 mm. Sketches of the BHE cross sections considered are illustrated in Figure 1. Water has been considered as working fluid. The water thermal properties at Tin have been taken from NIST [44]: density ρw = 995.03 kg/m 3 , dynamic viscosity μw = 0.76456 mPa s, specific heat capacity cpw = 4179.5 J/(kg K), thermal conductivity kw = 0.61869 W/(mK). The water velocity was considered vertical and uniform. A heat flux per unit area given by the product of the convection coefficient and the temperature difference between fluid and solid surface was applied at the fluid-solid interface. The Reynolds, Nusselt and Prandtl numbers, and the heat transfer coefficient h are given in Table 1. The Nusselt number was calculated through the Churchill correlation with uniform wall heat flux [45]. For each geometry, we considered three grout thermal conductivities, 1.0, 1.6 and 2.3 W/(mK), two flow rates, . V = 12 and 24 L/min, and ground thermal conductivity k g = 1.8 W/(mK) (typical value).
The following values were adopted for the thermal properties nearly uninfluential on ϕ: polyethylene thermal conductivity and volumetric heat capacity k p = 0.4 W/(mK) and (ρ c) p = 1.824 MJ/(m 3 K); grout and ground volumetric heat capacities (ρ c) gt = 1.600 MJ/(m 3 K) and (ρ c) g = 2.500 MJ/(m 3 K). We examined the cooling operation, with inlet fluid temperature 32 • C. Thus, 54 finite element simulations were performed to determine the correlations for ϕ. Additional simulations were performed to validate the simulation code, as well as to check the validity of the correlations for other values of D b and of k g and for other working conditions.
Water has been considered as working fluid. The water thermal properties at T in have been taken from NIST [44]: density ρ w = 995.03 kg/m 3 , dynamic viscosity µ w = 0.76456 mPa s, specific heat capacity c pw = 4179.5 J/(kg K), thermal conductivity k w = 0.61869 W/(mK). The water velocity was considered vertical and uniform. A heat flux per unit area given by the product of the convection coefficient and the temperature difference between fluid and solid surface was applied at the fluid-solid interface. The Reynolds, Nusselt and Prandtl numbers, and the heat transfer coefficient h are given in Table 1. The Nusselt number was calculated through the Churchill correlation with uniform wall heat flux [45].  The ground surrounding the borehole has been represented as a cylinder coaxial with the borehole, having radius 10 m and 10 m longer than the BHE. The initial temperature of the ground and of the BHE has been set equal to the undisturbed ground temperature, T g . The latter has been supposed equal to 14 • C for z = 10 m, with geothermal gradient 0.03 • C/m for z > 10 m. The following distribution of T g (z) has been assumed for z < 10 m: An adiabatic boundary condition has been imposed at the lateral and bottom ground surfaces and at the top of the borehole. The boundary condition T g (0) = 24 • C has been applied at the horizontal ground surface.
Numerical simulations for a working period of 100 h have been carried out by a 3D finite element model, through COMSOL Multiphysics. The working period selected is more than sufficient to reach a quasi-stationary heat transfer regime in the BHE and to obtain abundant data of ϕ in this regime.
The reduced vertical coordinate z = z/c has been introduced to shorten the computational domain along z [46]. Consequently, a reduced thermal conductivity along z, k z = k z /c 2 , has been employed for every material; in addition, a reduced vertical water velocity w = w/c has been assumed. A more detailed description of the method can be found in Refs. [43,46]. The reduction coefficient c has been taken equal to 5 for BHEs with length 50 m, equal to 10 for BHEs with length 100 m, and equal to 20 for BHEs with length 200 m.
In COMSOL Multiphysics, the time steps are non-uniform and are optimized by the software so that the solution matches the accuracy parameters imposed by the user. We have selected absolute tolerance 0.0001 and relative tolerance 0.001, instead of the default values 0.001 and 0.01, respectively.
For each geometry, the computational domain has been meshed with unstructured tetrahedral elements. Due to the rescaling coefficient, the selected mesh, illustrated in Figure 2, is independent of the BHE length, and has 1,369,572 elements for d = 84 mm, 1,444,394 elements for d = 94 mm, and 1,510,647 elements for d = 104 mm.  The ground surrounding the borehole has been represented as a cylinder coaxial with the borehole, having radius 10 m and 10 m longer than the BHE. The initial temperature of the ground and of the BHE has been set equal to the undisturbed ground temperature, Tg. The latter has been supposed equal to 14 °C for z = 10 m, with geothermal gradient 0.03 °C/m for z > 10 m. The following distribution of Tg(z) has been assumed for z < 10 m: An adiabatic boundary condition has been imposed at the lateral and bottom ground surfaces and at the top of the borehole. The boundary condition Tg(0) = 24 °C has been applied at the horizontal ground surface.
Numerical simulations for a working period of 100 h have been carried out by a 3D finite element model, through COMSOL Multiphysics. The working period selected is more than sufficient to reach a quasi-stationary heat transfer regime in the BHE and to obtain abundant data of ϕ in this regime.
The reduced vertical coordinate =  / z z c has been introduced to shorten the computational domain along z [46]. Consequently, a reduced thermal conductivity along z, =  2 z z k k c , has been employed for every material; in addition, a reduced vertical water velocity =  w w c has been assumed. A more detailed description of the method can be found in Refs. [43,46]. The reduction coefficient c has been taken equal to 5 for BHEs with length 50 m, equal to 10 for BHEs with length 100 m, and equal to 20 for BHEs with length 200 m.
In COMSOL Multiphysics, the time steps are non-uniform and are optimized by the software so that the solution matches the accuracy parameters imposed by the user. We have selected absolute tolerance 0.0001 and relative tolerance 0.001, instead of the default values 0.001 and 0.01, respectively.
For each geometry, the computational domain has been meshed with unstructured tetrahedral elements. Due to the rescaling coefficient, the selected mesh, illustrated in Figure 2, is independent of the BHE length, and has 1,369,572 elements for d = 84 mm, 1,444,394 elements for d = 94 mm, and 1,510,647 elements for d = 104 mm.  In order to ensure that the results are mesh independent, simulations have been carried out by employing meshes with 1,198,027, 1,311,663 and 1,444,394 tetrahedral elements, for L = 100 m, d = 94 mm, k gt = 1.6 W/(mK), k g = 1.8 W/(mK), and . V = 12 L/min. The values of T ave − T fm determined at τ = 2 h, 20 h and 100 h from the operation start have been compared, and the maximum percent deviation from the values obtained with the third mesh, employed in the final simulations, has been 0.127%.

Correlations for ϕ
The correlations for ϕ have been determined first for the quasi-stationary regime within the BHE, and then for the transient regime.

Quasi-Stationary Regime
The simulation results revealed that, after 2 h of operation with constant T in and . V, the quantity T ave − T fm becomes a homogeneous linear function of (T in − T out )/ . V that depends only on L, d, and k gt . Thus, a homogeneous linear regression of T ave − T fm as a function of (T in − T out )/ . V has been determined for each BHE geometry and for each value of k gt , by employing the simulation results for τ ≥ 2 h. The angular coefficient of each expression, divided by . V 0 = 12 L/min, gives a value of ϕ for the quasi-stationary regime, denoted by ϕ ∞ . The linear regressions obtained for L = 100 m and d = 94 mm are illustrated, as an example, in Figure 3. The values of ϕ ∞ are reported in Table 2.
Energies 2020, 13, x FOR PEER REVIEW 6 of 23 2 h, 20 h and 100 h from the operation start have been compared, and the maximum percent deviation from the values obtained with the third mesh, employed in the final simulations, has been 0.127%.

Correlations for φ
The correlations for ϕ have been determined first for the quasi-stationary regime within the BHE, and then for the transient regime.

Quasi-Stationary Regime
The simulation results revealed that, after 2 h of operation with constant Tin and V  , the quantity  Table 2. In order to determine a correlation for ϕ∞, the following dimensionless quantities have been employed: where L0 = 100 m, kgt0 = 1.6 W/(mK), and d0 = 94 mm are reference values of L, kgt, and d. By employing these parameters, we obtained the correlation:  In order to determine a correlation for ϕ ∞ , the following dimensionless quantities have been employed: where L 0 = 100 m, k gt0 = 1.6 W/(mK), and d 0 = 94 mm are reference values of L, k gt , and d. By employing these parameters, we obtained the correlation: Equation (11) can be employed to determine ϕ ∞ , and therefore T ave − T fm and T out − T fm in quasi-stationary regime, for every single U-tube borehole with 50 m ≤ L ≤ 200 m, 84 mm ≤ d ≤ 104 mm, 1.0 W/(mK) ≤ k gt ≤ 2.3 W/(mK). The mean square deviation of the values of ϕ ∞ given by Equation (11) from those reported in Table 2 is equal to 0.00237, i.e., to 2.51% of the mean value of ϕ ∞ .

Transient Regime
The simulation results revealed that the time dependence of ϕ during the two first working hours can be fitted by the equation: where a and b are dimensionless coefficients dependent on L, d, k gt and . V, and t * is the dimensionless time, evaluated as the working time τ divided by two hours. Moreover, the dependence on L and . V can be reduced to the dependence on the dimensionless parameter: Four values of V * were considered, namely 0.5, 1, 2 and 4. The coefficients a and b were determined by fitting values of ϕ/ϕ ∞ as a function of t * obtained numerically. Figure 4 illustrates the best fit for V * Energies 2020, 13, 828 8 of 23 = 1, d * = 1, k * = 1, that was obtained through two simulations: with L * = 1, As shown by the figure, the best-fit curve yields very accurate values of ϕ/ϕ ∞ , except for 0.125 < t * < 0.25, where it slightly overestimates the simulation results. The obtained values of a and b are listed in Table 3. As revealed by the table, a depends on V * , d * , and k * , while b depends only on V * . The correlation: fits the values of a given by Table 3 with a mean square deviation equal to 0.246. Table 3. Values of a and b to be used in Equation (12).

Figure 4. Best fit of the simulation results
The obtained values of a and b are listed in Table 3. As revealed by the table, a depends on V * , d * , and k * , while b depends only on V * . The correlation: fits the values of a given by Table 3 with a mean square deviation equal to 0.246. The relative discrepancy between Equation (14) and Table 3 can be considered as negligible for values of a higher than 10.
The coefficient b, which depends only on V * , can be expressed as: Equations (11), (12), (14) and (15) yield the time evolution of ϕ during the first two working hours, for every single U-tube BHE, in any working condition.
After the first hour of operation, the ratio ϕ/ϕ ∞ is very close to the asymptotic value 1. By integrating Equation (12) between t * = 0 and t * = 0.5 one obtains the mean value of ϕ during the first hour: Table 3. Values of a and b to be used in Equation (12).

Validation of the 3D Simulation Code
The validation of the 3D simulation code has been performed by comparison between the time-dependent values of T fm evaluated through this code and those calculated analytically through the BHE model proposed by Man et al. [35]. The BHE selected for the comparison has L = 100 m, d = 94 mm, k gt = 1.6 W/(mK), (ρ c) gt = 1.600 MJ/(m 3 K), volume flow rate 18 L/min, is placed in a ground with k g = 1.8 W/(mK) (ρ c) g = 2.500 MJ/(m 3 K), T g = 14 • C, and receives a time constant linear power of 60 W/m, i.e., to a total thermal power . Q = 6000 W. The initial temperature coincides with T g (z) in the whole domain. The properties of water are evaluated at 20 • C and are [44] ρ w = 998.21 kg/m 3 , µ w = 1.0016 mPa s, c pw = 4184.1 J/(kg K), k w = 0.59846 W/(mK). The convection coefficient, calculated by the Churchill correlation [45], is h = 1800.5 W/(m 2 K).
In the 3D simulation code, the condition of constant thermal power . Q = 6000 W has been implemented by imposing the constant value of T in − T out given by Equation (3), namely 4.789 • C. The mesh is that illustrated in Figure 2.
In the 2D axisymmetric BHE model proposed by Man et al. [35], the BHE is represented by a cylindrical surface with radius r 0 that releases a constant and uniform power per unit area corresponding to the total thermal power received by the BHE. The generating surface represents the fluid, and the borehole has the same properties as the ground. The analytical solution for the temperature field is determined by the Green's function method.
The mean temperature of the surface with radius r 0 , that will be denoted by T fm , is given by [35]: where α g = k g /(ρc) g is the thermal diffusivity of the ground, equal to 0.72 × 10 −6 m 2 /s in the case considered.
The value of r 0 to be employed in the model has been determined by imposing that the thermal resistance of the cylindrical layer between r 0 and the BHE radius is equal to the BHE thermal resistance. The latter has been determined through a stationary 2D finite element simulation of a borehole cross section that includes a ground layer having radius 2 m and an isothermal external surface. The temperature difference between the fluid and the external ground surface has been set equal to 25 • C for one pipe and to 20 • C for the other pipe. The convection coefficient is h = 1800.5 W/(m 2 K). We adopted a very fine mesh, composed of 139,008 triangular elements. A particular of the mesh employed is illustrated in Figure 5. The result is R b = 0.09863 mK/W and, as a consequence, r 0 = 2.491 cm. The value of r0 to be employed in the model has been determined by imposing that the thermal resistance of the cylindrical layer between r0 and the BHE radius is equal to the BHE thermal resistance. The latter has been determined through a stationary 2D finite element simulation of a borehole cross section that includes a ground layer having radius 2 m and an isothermal external surface. The temperature difference between the fluid and the external ground surface has been set equal to 25 °C for one pipe and to 20 °C for the other pipe. The convection coefficient is h = 1800.5 W/(m 2 K). We adopted a very fine mesh, composed of 139,008 triangular elements. A particular of the mesh employed is illustrated in Figure 5. The result is Rb = 0.09863 mK/W and, as a consequence, r0 = 2.491 cm. Figure 5. Particular of the mesh adopted to determine Rb, for the validation of the 3D simulation code. The boundary between the polyethylene pipes and the grout is hidden by the triangular elements.
The diagrams of Tfm versus the decimal logarithm of time in hours obtained by the 3D finite element simulation and by the numerical integration of Equation (17) are compared in Figure 6, in the time range between 10 −2 h and 10 2 h. The mean square deviation between the results of the finite element simulation and those obtained by the model by Man et al. [35] is 0.26 °C. The small discrepancies that occur from 10 −2 to 10 −1 h are very probably due to the non-perfect accuracy of the analytical model, that does not consider the exact total value and distribution of the BHE heat capacity. On the contrary, those occurring from 10 to 10 2 h are probably due to numerical errors in the 3D simulation, that increase with time. The diagrams of T fm versus the decimal logarithm of time in hours obtained by the 3D finite element simulation and by the numerical integration of Equation (17) are compared in Figure 6, in the time range between 10 −2 h and 10 2 h. The mean square deviation between the results of the finite element simulation and those obtained by the model by Man et al. [35] is 0.26 • C. The small discrepancies that occur from 10 −2 to 10 −1 h are very probably due to the non-perfect accuracy of the analytical model, that does not consider the exact total value and distribution of the BHE heat capacity. On the contrary, those occurring from 10 to 10 2 h are probably due to numerical errors in the 3D simulation, that increase with time.
the time range between 10 −2 h and 10 2 h. The mean square deviation between the results of the finite element simulation and those obtained by the model by Man et al. [35] is 0.26 °C. The small discrepancies that occur from 10 −2 to 10 −1 h are very probably due to the non-perfect accuracy of the analytical model, that does not consider the exact total value and distribution of the BHE heat capacity. On the contrary, those occurring from 10 to 10 2 h are probably due to numerical errors in the 3D simulation, that increase with time.

Validity of the Correlations for Other BHE Diameters, Thermal Conductivities of the Ground, Working Conditions
Although the correlations for ϕ reported in Section 3 were obtained by assuming D b = 152 mm, k g = 1.8 W/(mK), and summer operation with constant T in , they hold also for other BHE diameters, ground thermal conductivities, and working conditions.
The applicability of the values of ϕ ∞ reported in Table 2 to other borehole diameters is shown in Figure 7, that refers to BHEs with L = 100 m, d = 84 mm, k gt = 1.0 W/(mK), k g = 1.8 W/(mK), . V = 12 L/min, and BHE diameters 134 mm and 170 mm. Clearly, the value of D b has no effect on ϕ if k gt = k g . Therefore, the value k gt = 1.0 W/(mK) has been selected, to analyze a critical condition, with a high ratio k g /k gt . Higher values of this ratio should be avoided. The diagrams of T ave − T fm versus (T in − T out )/ . V obtained by the 3D simulations with these BHE diameters are compared with the line obtained by employing the value of ϕ ∞ reported in Table 2, namely ϕ ∞ = 0.07. The comparison shows the applicability of the correlation to other BHE diameters.
The applicability of the values of ϕ ∞ reported in Table 2 to other ground thermal conductivities is shown in Figure 8 Table 2, namely ϕ ∞ = 0.0847. The results reveal that the correlation can be applied to other values of the ground conductivity.
The applicability of the values of ϕ ∞ reported in Table 2 to other working conditions is shown in Figure 9, that refers to thermal response tests (TRTs) performed on a borehole with L = 100 m, d = 94 mm, k gt = 1.6 W/(mK), k g = 1. V is obtained for each TRT, in the quasi-stationary regime. In Figure 9, the values of T ave − T fm as a function of (T in − T out )/ . V obtained for the TRTs described above are compared with the line obtained by employing the value of ϕ ∞ reported in Table 2. The points that represent the results for the TRTs lay on the correlation line.
The applicability of the values of ϕ∞ reported in Table 2 to other ground thermal conductivities is shown in Figure 8, that refers to BHEs with L = 100 m, d = 94 mm, kgt = 1.6 W/(mK), V  = 12 L/min,  Table 2, namely ϕ∞ = 0.0847. The results reveal that the correlation can be applied to other values of the ground conductivity.   Table 2.  Table 2.
The applicability of the values of ϕ∞ reported in Table 2 to other working conditions is shown in  Table 2. The points that represent the results for the TRTs lay on the correlation line. V obtained by 3D simulations, for k g = 1.4 W/(mK) and k g = 2.2 W/(mK), and by the value of ϕ ∞ reported in Table 2. T T V is obtained for each TRT, in the quasistationary regime. In Figure 9, the values of Tave − Tfm as a function of ( ) −  in out T T V obtained for the TRTs described above are compared with the line obtained by employing the value of ϕ∞ reported in Table 2. The points that represent the results for the TRTs lay on the correlation line. V obtained for TRTs, compared with the correlation given by the value of ϕ ∞ reported in Table 2.

Evaluation of the Time Evolution of T out by Means of the Correlations for ϕ
In this section, we show how the correlations for ϕ determined in Section 3 can be used to evaluate accurately T out by means of simulation codes that yield T fm , first in the case of constant heat load and flow rate, then in the case of time dependent heat load and constant flow rate. The accuracy obtainable in working conditions with time dependent flow rate is not analyzed in this paper.

Constant Heat Load
In this example, we determine T fm by the analytical BHE model by Man et al. [35] and T out by the correlations for ϕ, and we point out the errors that occur in the short time if one employs simpler methods, such as the finite line-source (FLS) model combined with either R b or R beff , or the model by Man et al. [35] without the correlations for ϕ.  [45]. The BHE thermal resistance, determined by a 2D steady state simulation by the method explained in Section 4, is R b = 0.09965 mK/W. The effective BHE thermal resistance has been determined by Hellström's analytical expression [47]: where: and R a is the internal resistance between the tubes. The latter has been evaluated through the expression [48]: where T f1 and T f2 are the bulk temperatures in the tubes, and q l1 and q l2 are the heat fluxes per unit length from the tubes. The results obtained are R a = 0.47555 mK/W, R b,eff = 0.10950 mK/W.
In the model of Ref. [35], the radius of the generating surface such that the thermal resistance of the cylindrical layer between that surface and the BHE radius equals R b is r 0 = 2.462 cm.
The simplest method to obtain the time evolution of T out is determining T sm by the FLS model and T fm by Equation (1), then adopting the approximation: T fm ≈ T ave (21) and evaluating T out by the relation: An improvement of this method is determining T sm by the FLS model, T ave by the definition of effective BHE thermal resistance: and then T out by Equation (22). This improvement allows obtaining correct values of T out in the long term.
The time evolution of T out obtained by applying the FLS model and R b and that obtained by applying the FLS model and R b,eff are compared in Figure 10 with that obtained by a 3D finite element simulation performed by the method described in Section 2 and the grid illustrated in Figure 2. The time interval considered is from 10 −2 h to 10 2 h, in the logarithmic scale. The time evolution of Tout obtained by applying the FLS model and Rb and that obtained by applying the FLS model and Rb,eff are compared in Figure 10 with that obtained by a 3D finite element simulation performed by the method described in Section 2 and the grid illustrated in Figure 2 The figure shows that the methods based on the FLS model yield an important overestimation of Tout during the first hour. In fact, while the true value of Tout remains practically equal to 14 °C during the first ten minutes, Tout jumps immediately from 14 to 16 °C according to the FLS_Rb method, and from 14 to 16.48 °C according to the FLS_Rb,eff method, that is more accurate in the long term. The overestimation of Tout by the FLS_Rb,eff method starts decreasing after 10 min, but remains relevant for about one hour (0.76 °C for τ = 1 h).
A much more accurate evaluation of the time evolution of Tout can be obtained by determining the time evolution of Tfm by means of a full-time scale BHE model that takes into account the heat capacity of the borehole elements, and calculating Tout through the coefficient ϕ and Equation (6). The time evolution of Tout yielded by the model of Ref. [35] and the coefficient ϕ is compared with that obtained by the 3D finite element simulation described above in Figure 11. The time dependent value of ϕ has been determined through Equation (12), with the value of ϕ∞ given in Table 2 (0.0847) and the values of a and b given in Table 3 (4.1 and 14, respectively). The figure shows that the combined use of the model by Man et al. [35] and of the coefficient ϕ yields results in fair agreement with those of the 3D numerical simulation, even for very short times.  Figure 11 illustrates also the results obtained by employing the model of Ref. [35] coupled with the approximation Tfm ≈ Tave and Equation (22), and shows that the approximate method underestimates Tout in the short time. The time evolution of T sm by the FLS model has been obtained through the numerical integration of the expression determined by Bandos et al. [15]: The figure shows that the methods based on the FLS model yield an important overestimation of T out during the first hour. In fact, while the true value of T out remains practically equal to 14 • C during the first ten minutes, T out jumps immediately from 14 to 16 • C according to the FLS_R b method, and from 14 to 16.48 • C according to the FLS_R b,eff method, that is more accurate in the long term. The overestimation of T out by the FLS_R b,eff method starts decreasing after 10 min, but remains relevant for about one hour (0.76 • C for τ = 1 h).
A much more accurate evaluation of the time evolution of T out can be obtained by determining the time evolution of T fm by means of a full-time scale BHE model that takes into account the heat capacity of the borehole elements, and calculating T out through the coefficient ϕ and Equation (6). The time evolution of T out yielded by the model of Ref. [35] and the coefficient ϕ is compared with that obtained by the 3D finite element simulation described above in Figure 11. The time dependent value of ϕ has been determined through Equation (12), with the value of ϕ ∞ given in Table 2 (0.0847) and the values of a and b given in Table 3 (4.1 and 14, respectively). The figure shows that the combined use of the model by Man et al. [35] and of the coefficient ϕ yields results in fair agreement with those of the 3D numerical simulation, even for very short times.

Time Dependent Heat Load
We illustrate now the use of our correlation for ϕ to determine Tout in the case of a heat load with sharp step changes. We consider the hourly demand of thermal energy of an uninsulated office building located in Bologna (Italy), during the first two days of the heating period, namely October 15 and 16, of a typical meteorological year. The building has two floors and a heated floor area of 1255.7 m 2 . Heating is turned on at 5 a.m. and is turned off at 7 p.m. The hourly values of the thermal energy required by the building, including the distribution and emission heat losses, have been determined by a dynamic simulation performed through TRNSYS 17. In our example, we assume that the building is heated by a heat pump with COP = 4, coupled with 12 BHEs, each 100 m long, and determine the hourly heat loads on each BHE multiplying by 0.75 and dividing by 12 the hourly values of the total thermal energy required by the building. A dynamic simulation of the buildingplant system with values of the heat pump COP determined at each hour is beyond the aims of this paper. The absolute values of the negative hourly heat load on each BHE are illustrated in Figure 12.  Figure 11 illustrates also the results obtained by employing the model of Ref. [35] coupled with the approximation T fm ≈ T ave and Equation (22), and shows that the approximate method underestimates T out in the short time.

Time Dependent Heat Load
We illustrate now the use of our correlation for ϕ to determine T out in the case of a heat load with sharp step changes. We consider the hourly demand of thermal energy of an uninsulated office building located in Bologna (Italy), during the first two days of the heating period, namely October 15 and 16, of a typical meteorological year. The building has two floors and a heated floor area of 1255.7 m 2 . Heating is turned on at 5 a.m. and is turned off at 7 p.m. The hourly values of the thermal energy required by the building, including the distribution and emission heat losses, have been determined by a dynamic simulation performed through TRNSYS 17. In our example, we assume that the building is heated by a heat pump with COP = 4, coupled with 12 BHEs, each 100 m long, and determine the hourly heat loads on each BHE multiplying by 0.75 and dividing by 12 the hourly values of the total thermal energy required by the building. A dynamic simulation of the building-plant system with The time evolutions of Tfm and of Tout are determined by a 3D finite element simulation of the BHE, performed by the same simulation code and the same mesh as in Section 6.1. Then, the time evolution of Tout is calculated from that of Tfm by means of our correlations for ϕ and Equation (6), with the following method. The time interval Δt needed by the fluid to go from the inlet to the outlet is evaluated by assuming plug flow, i.e., a uniform velocity in each cross section: 14 min for the case considered. For each hour, the dimensionless time t * of Equation (12) starts from zero at the beginning of the hour, and the application of the correlations for ϕ starts after a time interval Δt from the beginning of the hour. During the time interval Δt of the first hour, Tout is kept constant with the same value as in the initial instant, as is physically consistent and confirmed by the 3D simulation. During the time interval Δt of any other hour, Tout is considered as linearly varying with time from the last time instant of the preceding hour to the first time instant after Δt of that hour. A comparison between the time evolution of Tout yielded by the 3D simulation and that obtained through our correlation for ϕ is illustrated in Figure 13, and a particular for the first five hours is illustrated in Figure 14. The figures report also the time evolution of Tout that is obtained by employing the same method to manage the initial time interval Δt of each hour, and by applying the customary approximation Tfm = Tave instead of the correlations for ϕ. The time evolutions of T fm and of T out are determined by a 3D finite element simulation of the BHE, performed by the same simulation code and the same mesh as in Section 6.1. Then, the time evolution of T out is calculated from that of T fm by means of our correlations for ϕ and Equation (6), with the following method. The time interval ∆t needed by the fluid to go from the inlet to the outlet is evaluated by assuming plug flow, i.e., a uniform velocity in each cross section: 14 min for the case considered. For each hour, the dimensionless time t * of Equation (12) starts from zero at the beginning of the hour, and the application of the correlations for ϕ starts after a time interval ∆t from the beginning of the hour. During the time interval ∆t of the first hour, T out is kept constant with the same value as in the initial instant, as is physically consistent and confirmed by the 3D simulation. During the time interval ∆t of any other hour, T out is considered as linearly varying with time from the last time instant of the preceding hour to the first time instant after ∆t of that hour. A comparison between the time evolution of T out yielded by the 3D simulation and that obtained through our correlation for ϕ is illustrated in Figure 13, and a particular for the first five hours is illustrated in Figure 14. The figures report also the time evolution of T out that is obtained by employing the same method to manage the initial time interval ∆t of each hour, and by applying the customary approximation T fm = T ave instead of the correlations for ϕ.  The figures show that the time evolution of Tout obtained by applying the correlations for ϕ is much more accurate, especially during the most critical hours. The mean square deviation from the results of the 3D simulation, during the first three hours, is 0.1 °C for the method employing our correlations for ϕ and 0.6 °C for that employing the approximation Tfm = Tave.  The figures show that the time evolution of Tout obtained by applying the correlations for ϕ is much more accurate, especially during the most critical hours. The mean square deviation from the results of the 3D simulation, during the first three hours, is 0.1 °C for the method employing our correlations for ϕ and 0.6 °C for that employing the approximation Tfm = Tave. The figures show that the time evolution of T out obtained by applying the correlations for ϕ is much more accurate, especially during the most critical hours. The mean square deviation from the results of the 3D simulation, during the first three hours, is 0.1 • C for the method employing our correlations for ϕ and 0.6 • C for that employing the approximation T fm = T ave .

Calculation of R b,eff Through ϕ ∞
In this section we show that, in quasi-stationary working conditions, the effective BHE thermal resistance, R b,eff , can be easily calculated by means of the dimensionless coefficient ϕ ∞ . Then, we compare the values of R b,eff evaluated through our general correlation for ϕ ∞ , given by Equation (11), with those determined by 3D finite element computations and by Hellström's analytical solution.
One can rewrite Equation (23) as: In quasi-stationary working conditions, the energy balance of the fluid yields: By substituting Equation (26) in Equation (25), one gets: In quasi-stationary conditions, Equation (5) yields: By substituting Equation (28) in Equation (27), one obtains: Equations (29) and (11) allow a very simple calculation of the difference between R b,eff and R b , for single U-tube BHEs. Equation (29) can be employed also in combination with Equation (27) of Ref. [43], for double U-tube BHEs. The BHE thermal resistance R b can be calculated either by employing a suitable approximate expression, or, more accurately, through a 2D numerical simulation of a borehole cross section.
A comparison between the values of R b,eff obtained by Equations (29) and (11), denoted by R b,eff,ϕ , those evaluated through 3D numerical computations, denoted by R b,eff,3D , and those obtained by Equations (18) and (19), denoted by R b,eff,H , is illustrated in Table 4 Table 1. The values of R a and of R b to be employed in Equations (18) and (19) have been determined by performing stationary 2D numerical simulations of a borehole cross section with the method described in Section 4 and evaluating R a with Equation (20).
The table shows that the values of R b,eff evaluated through our correlation for ϕ ∞ are in fair agreement both with those computed directly by 3D numerical computations and with those yielded by Hellström's analytical method. The highest percent deviation from R b,eff,3D is −1.46% for R b,eff,ϕ and −1.81% for R b,eff,H ; the mean square deviation from R b,eff,3D is 0.00093 mK/W for R b,eff,ϕ and 0.00119 mK/W for R b,eff,H .
The values of R b,eff,ϕ reported in Table 4 have been determined by calculating ϕ ∞ through Equation (11). If one employs the values of ϕ ∞ reported in Table 2, the highest percent deviation from R b,eff,3D reduces to −1.05% and the mean square deviation reduces to 0.00086 mK/W.

Conclusions
By means of 3D finite element simulations of single U-tube borehole heat exchangers (BHEs), we have obtained simple expressions of a dimensionless coefficient ϕ that yields two temperature differences: that between the average of inlet and outlet fluid temperature, T ave , and the mean fluid temperature, T fm ; that between the outlet temperature, T out , and T fm . The obtained expressions of ϕ hold for single U-tube BHEs with any diameter between 130 and 170 mm, any shank spacing compatible with the diameter, and any BHE length between 50 m and 200 m, both in quasi-stationary and in transient regime. The 3D simulation code has been validated by comparing the time evolution of T fm obtained by the code with that obtained by applying an analytical BHE model.
We have shown that the obtained expressions of ϕ allow an accurate evaluation of the time evolution of T out through a BHE model that yields the time evolution of T fm , for working conditions with constant flow rate and either constant or time dependent heat load. Moreover, we have illustrated the errors that occur in the short time if one applies the finite line-source model coupled with the effective BHE thermal resistance, and those that occur if one applies a BHE model that yields T fm coupled with the approximation T fm = T ave . Thus, the obtained expressions of ϕ are a useful complement of the BHE models that yield T fm . Finally, we have shown that the obtained quasi-stationary values of ϕ can be used for a simple evaluation of the effective BHE thermal resistance.
Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.