Temperature response factors at different boundary conditions for modelling the single borehole heat exchanger

Applied Thermal Engineering 103 (2016) 934–944

Contents lists available at ScienceDirect

Applied Thermal Engineering journal homepage: www.elsevier.com/locate/apthermeng

Research Paper

Temperature response factors at different boundary conditions for modelling the single borehole heat exchanger Antonella Priarone ⇑, Marco Fossa DIME, University of Genova, via Opera Pia 15a, 16145 Genova, Italy

h i g h l i g h t s
Numerical model to calculate the Temperature Response Factors for a single BHE.
Different BHE wall boundary conditions analyzed.
Model able to consider the BHE position with respect to the ground surface.
Obtained results provided in a tabulated form.
Data useful for spatial and temporal superposition for g-function generation.

a r t i c l e

i n f o

Article history: Received 2 March 2016 Revised 7 April 2016 Accepted 8 April 2016 Available online 9 April 2016 Keywords: Ground coupled heat pumps Borehole heat exchanger design Temperature response factors g-functions

a b s t r a c t Design and simulation of borehole heat exchangers rely on the solution of the transient conduction equation. The typical approach for predicting the ground temperature variations in the short and long term is to recursively apply basic thermal response factors available as analytical functions or as preestimated tabulated values. In this paper a review of the existing response factor models for borehole heat exchangers (BHE) analysis is presented and a numerical model, built in Comsol environment is employed for calculating the temperature distribution in time and space around a single, finite length, vertical cylindrical heat source also taking into account its position with reference to the ground surface (effects of the adiabatic length or ‘‘buried depth” D). The temperature values are recast as dimensionless response factors in order to compare them with analytical solutions where available. Furthermore new temperature response factors suitable for describing the single Finite Cylindrical Source (FCS) under different operating modes (i.e. boundary conditions) are generated. Boundary conditions include imposed heat transfer rate, imposed temperature and a combination of both conditions, where spatially uniform temperature at the BHE interface is attained while also keeping constant the applied heat transfer rate. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction Ground coupled heat pumps (GCHP) are probably the most energy efficient solution for building space conditioning. This technology has encountered a wide diffusion in several cold climate countries, either in Northern Europe or North America where traditional air-to-air heat pumps cannot be efficiently employed for heating purposes at very low ambient (outdoor) air temperatures. This technology can effectively concur to the attainment of the EU energy and climate targets. To promote a widespread exploitation of these resources on field and a large scale introduction of GCHP systems in both residential and institutional/commercial buildings, it is mandatory to define either technical guidelines for ⇑ Corresponding author. E-mail address: [email protected] (A. Priarone). http://dx.doi.org/10.1016/j.applthermaleng.2016.04.038 1359-4311/Ó 2016 Elsevier Ltd. All rights reserved.

the correct sizing of ground heat exchanger fields or a set of regulatory and economic actions [1]. The most common solution for extracting or injecting heat from and to the ground is to bury closed loop heat exchangers either arranged near the ground surface or disposed vertically down to hundreds of meters (BHEs). If reliable drilling equipment is available, vertical heat exchangers are usually preferred to near horizontal or trench siblings due to the reduced requirement of land surface and to their capability of taking advantage of the stable and favorable temperatures of the deep soil. GCHP can cover a wide range of energy demand situations, from small residences to large commercial buildings. High efficiencies in such installations are related to the correct design of the borehole field in order to have suitable return temperatures of the heat carrier fluid for achieving high seasonal performance factors. Such a task can be accomplished through the correct modelling of the ground/

A. Priarone, M. Fossa / Applied Thermal Engineering 103 (2016) 934–944

935

Nomenclature B c D erf erfc E1 F Fo g G H J0, J1 k ~ n p Q_ r T Y0, Y1 z

radial distance from multiple sources (m) specific heat (J/(kg K)) adiabatic depth (m) error function (–) complementary error function (–) exponential integral in ILS model (–) Infinite Cylindrical Source (ICS) temperature transferfunction for imposed temperature (–) Fourier number (–) multi BHE temperature transfer function (–) ICS temperature transfer-function for imposed heat flux (–) active depth of the BHE (m) Bessel functions of the first kind, of order 0 and 1 (–) thermal conductivity (W/(m K)) inward normal unit vector (–) ratio r/rb (–) heat transfer rate (W) radial coordinate (m) temperature (K) Bessel functions of the second kind, of order 0 and 1 (–) vertical coordinate (m)

borefield system and suitable models are required to simulate the short to long time temperature evolution of the ground volume associated with the BHE field. From a thermal point of view, this volume can be very often considered as a purely conduction medium since ground water circulation (if present) is often confined in a thin layer compared to the vertical extension of the BHE volume. When meaningful underwater circulation is present a common approach is to include convection in the conduction equation (see for example [2] but often a pure conduction modelling, in terms of moving line sources, demonstrated to be able to efficiently describe the ground thermal behavior in presence of groundwater advection [3]. The borefield design goal is the definition of the best BHE geometry (BHE arrangement, their number and spacing) and the minimum overall length of vertical pipes with respect to the land availability and drilling equipment. The constraints of the problem and its starting information are the building heat loads in time, the ground thermal properties and the target seasonal heat pump performance. In order to succeed in this design task, a number of solutions of the transient heat conduction equation have been proposed in order predict the ground temperature in time and space for given geometry and boundary conditions. Carslaw and Jaeger [4] and Ingersoll et al. [5] first provided engineering models for evaluating the temperature evolution in the ground when a single heat source is active. In both these early studies the borehole is modelled as an infinite length heat source. The most popular solutions for such a problem are the so called Infinite Line Source (ILS, Lord Kelvin, and later Ingersoll et al.) and Infinite Cylindrical Source (ICS, Carslaw and Jaeger). Both solutions allow the temperature distribution in the ground to be evaluated in terms of a temperature response factor (TRF) which is a function of the radius based Fourier number. ILS and ICS models were both solved for the constant heat flux boundary condition and the ICS for the imposed wall temperature too. Imposed heat rate cases result in similar trends for ILS and ICS models except for the early part of the temperature evolution. These two solutions can be proved

Greek letters thermal diffusivity (m2/s) density (kg/m3) time (s) C temperature transfer function (–)

a q s

Subscripts ave average b of the borehole, at borehole radius gr of the ground medium, of the ground domain H based on BHE active depth p r/rb Q imposed heat flux r based on radius sc superconductive material T imposed temperature 1 far field and initial condition Superscripts ⁄ dimensionless  average along the depth H 0 per unit length

to be in absolute agreement after a given dimensionless time is elapsed as also discussed in recent papers [6,7]. The ICS and ILS solutions at imposed heat rate encountered a great success in many engineering models, from short time analyses (TRT experiments, e.g. [8–11]) to long term simulations [12]. The ICS solution at imposed temperature, also known as the ‘‘F” function is rarely employed probably due to the fact that in practice this boundary condition is difficult to realize in real experiments or field equipment. Thanks to the work of the Lund research group (e.g. [13], the TRF approach was extended to the description of finite length (linear) heat sources (FLS). Lamarche and Beauchamp [14], after Zeng et al. [15], developed a new solution for the Finite Line Source model (FLS) in terms of semi analytical expressions to be numerically integrated. Bandos et al. [16] proposed new analytical and explicit expressions for calculating the ground temperature field induced by a finite line source. Fossa [17] and Fossa and Rolando [18] refined the Bandos expressions and proposed approximate fully analytical solutions for the FLS problem suitable for spatial and temporal superposition in very reduced computation times. Claesson and Javed [19] rearranged the mathematical development of the (semi analytical) FLS solution in a way able to provide new expression which greatly reduces the computation time of the response factor values. One of the strength of the TRF approach is the possibility to exploit the linear properties of the conduction Fourier equation for superposing the base solution in space and hence to obtain new TRFs able to describe the ground response to a system of heat sources. These multiple source TRFs are usually referred as ‘‘gfunctions” after the work of Eskilson [13]. The approach of the Lund group was to numerically calculate the temperature field from a single (finite) heat source either at imposed heat flux or imposed temperature. The code in charge to manage such a problem, discretizing the Fourier equation according to a finite difference scheme, was named SBM [20]; its g-functions are embedded in commercial codes for BHE field design as EED and GLHEPRO [21,22]. From the point of view of numerical approach, Zanchini and Lazzari [23] proposed a new method that considers also the inter-

936

A. Priarone, M. Fossa / Applied Thermal Engineering 103 (2016) 934–944

nal structure of the BHE and provides new g-functions that yield the dimensionless temperature at the interface tubes-grout. Following the Eskilson [13] contribution, a number of authors applied the response factor approach and the superposition techniques to generate semi-analytical g-functions (i.e. based on semi-analytical FLS solutions), thus avoiding to perform time consuming and challenging numerical solutions of the discretized Fourier equation [24–27]. In spite of the general agreement among the FLS generated g-functions and the numerically derived SBM ones, discrepancies were noticed by all the above mentioned authors. These differences, noticeable when the BHE number is high, say for dozens of heat exchanges, have been ascribed to the different boundary conditions which pertain to the FLS solutions and the SBM ones. This problem was recently discussed by Cimmino et al. [28], who proposed an algorithm for calculating gfunctions referring to constant temperature BHEs through a modulation of the heat transfer rate at each individual BHE. In this paper the ILS, ICS and FLS models are described and a numerical model in Comsol environment [29] is further developed [30] to solve the heat conduction problem according to three different boundary conditions, namely the imposed heat flux condition, the imposed temperature condition and a combination of the previous two, when a spatially constant temperature is attained while still keeping a constant in time heat transfer rate, as recently also done by Monzo and Acuna [31] and by Monzo et al. [32]. Furthermore the Comsol model has been modified to take into account the presence of an adiabatic top part of the BHE, as possible in the SBM code and also made possible by applying the FLS formulas in the Claesson and Javed [19] version. The results are here presented with reference to the single heat source, while a next paper will be devoted to a complete analysis of the superposed numerical solutions (e.g. the corresponding g-functions) when different boundary conditions apply, thus contributing to a more complete view to the problem of the numerical and analytical generation of single and multiple source temperature response factors.

2. The single heat source theory applied to the GCHP problem When the underground water circulation can be neglected, the heat transfer between the ground and a vertical borehole heat exchanger can be described in terms of the three-dimensional time-dependent conduction equation. Depending on the assumptions related to the BHE geometry, the problem can be managed according to either a one-dimensional (in the radial direction) scheme or a two-dimensional (radial and axial) one. The simplest model is the one known as the Infinite Line Source (ILS), which was clearly discussed by Ingersoll et al. [5] as a possible mathematical solution for geothermal analyses. The ILS expression is based on a heat flux (per unit length Q_ 0 ) boundary condition:

Tðr; sÞ  T gr;1 ¼

sions referred to the active length H of the finite heat source and including the multiple evaluation of the erf function: T av e ðr; sÞ  T gr;1 ¼

"Z # 1 Q_ 0 2 2 Yðz; D=H  zÞ exp½ðr=HÞ z   dz z2 4pkgr pffiffiffiffiffiffiffiffiffiffiffi 1=4FoH

Yðx; yÞ ¼ 2ierf ðxÞ þ 2ierf ðx þ 2yÞ  ierf ð2x þ 2yÞ  ierf ð2yÞ Z U 2 1 ierf ðUÞ ¼ erf ðv Þdv ¼ U erf ðUÞ  pffiffiffiffi ð1  eU Þ

p

o

ð2Þ ð3Þ ð4Þ

The equation set (2)–(4) also includes the ‘‘buried depth” D, which represents the distance from the ground surface to the beginning of the line source of length H. It can also be considered as the adiabatic (top) part of a BHE having an overall length equal to (D + H). This additional parameter was already conceived by Eskilson [13] and made available as calculating option in the SBM code. Worth noticing, the equation set proposed by Lamarche and Beauchamp [14] and Eqs. (2)–(4) yield the same results when D is set to zero, but, as discussed by Fossa and Rolando [18], the most recent equation set allows faster calculations. Another group of solutions is the one belonging to the infinite cylindrical source family. This model was first proposed by Carslaw and Jaeger [4] and later described by Ingersoll et al. [5] who also provided tabulated values of the related solutions. The ICS case was analytically solved at imposed heat transfer rate and also at imposed temperature. After Carlslaw and Jaeger, the corresponding solutions of the above problems are referred as the ‘‘G” and ‘‘F” functions, respectively. They can be expressed according the formulas:   r Q_ 0 Tðr; sÞ  T gr;1 ¼ G Forb ; p ¼ rb kgr Z 1 b2 Forb 1 e 1 1 ¼ 2 ½J ðpbÞY 1 ðbÞ  J1 ðbÞY 0 ðpbÞ 2 db ð5Þ p 0 J21 ðbÞ þ Y 21 ðbÞ 0 b Z 2 8 1 eb Forb 1 db ð6Þ Q_ 0 ðsÞ ¼ kgr ðTðr b Þ  T gr;1 ÞFðForb Þ ¼ 2 p 0 J0 ðbÞ þ Y 20 ðbÞ b where J0, J1, Y0, Y1 are Bessel functions of the zeroth and first order, respectively. All the above expressions can be summarized introducing a generic temperature response factor C:

C ¼ ðT  T gr;1 Þ

 Q_ 0 2pkgr

ð7Þ

which can be related to either a constant temperature condition or to a constant heat transfer rate one. In the former case, the heat transfer rate (per unit length) is the average along the depth while in the latter case is in turn the temperature to be averaged along the heat source. 3. BHE and ground volume modelling 3.1. Equation set and non-dimensionalization

Z 1 b e Q_ 0 Q_ 0 E1 ð1=4For Þ db ¼ 4pkgr 4Fo1 b 4pkgr

ð1Þ

r

Here kgr is the ground thermal conductivity, Tgr,1 is the undisturbed (initial) ground temperature, For is the radius based Fourier number. The E1 function is referred as the exponential integral, which has a very useful expression as expansion in series. To describe the effects of a finite length line source Eskilson [13] and later Zeng et al. [15] introduced the Finite Line Source (FLS) approach, then refined by Lamarche and Beauchamp [14], who proposed a particular expression of the FLS solution involving the complementary error function (erfc). More recently Claesson and Javed [19] reformulated the FLS theory according to new expres-

A numerical model has been built in COMSOL Multiphysics environment in order to investigate the behavior of single borehole heat exchangers as described by Finite Cylindrical heat Sources (FCS) having different boundary conditions. The detailed description of the model, especially regarding model validation against the mesh characteristics and computational domain extension, has been largely discussed in a recent paper by the same Authors [30] and the reader is addressed to this work for further information on the present model. The computational domain refers to the ground where constant thermal properties (qgr, cgr and kgr, density, specific heat and thermal conductivity, respectively) are assumed. The BHE is represented by an inner surface of the axial–symmetrical domain,

937

A. Priarone, M. Fossa / Applied Thermal Engineering 103 (2016) 934–944

located at the radius r = rb. Pure heat conduction, in transient or steady state mode is considered. The transient general equation is hence:

@T @ 2 T 1 @T @ 2 T qgr cgr ¼ kgr r2 T ¼ kgr þ þ @s @r 2 r @r @z2

rgr rb

! ð8Þ

D

The ground volume is a large cylinder of radius rgr and depth Hgr. It is intended to be confined at its top by the external environment where an isothermal and constant in time condition is set. This temperature is also the undisturbed ground temperature Tgr,1, which represents either the initial condition of the problem or the boundary condition at the far end of the domain. If the domain is large enough with respect to the transient period under investigation, the far end boundary condition can be replaced by an adiabatic condition, as also discussed in Priarone and Fossa [30]. The BHE surface has a thermally active depth equal to H; the heat source can be also located a little bit far from top surface by a vertical displacement D. The possibility of this displacement was considered by the Lund research group that named it the ‘‘buried depth”. Probably a better definition would be the ‘‘adiabatic depth”, so indicating the BHE top part which is not active from a thermal point of view (e.g. in Swedish applications, where the borehole top part is not filled by water or it is deliberately insulated). On the BHE surface (i.e. the inner surface at r = rb) three different boundary conditions have been considered. The results pertaining two of them have been partially presented in the above mentioned recent paper by present Authors. The three boundary conditions at the BHE surface are: constant heat transfer rate per unit length (Eq. (9a)), constant temperature (Eq. (9b)), constant heat transfer rate at the inner surface of a superconductive core virtually inserted into the BHE hole (Eq. (9c)). The three case studies are depicted in Fig. 1. The third boundary condition has been realized by inserting an additional hollow cylinder, concentric with the BHE with external and internal radius equal to rb and rsc, respectively. The inner diameter rsc is set to 0.1rb. The inner core is made by a ideal thermally super-conductive material characterized by the following values

H

Q_ 0 ðD < z < D þ H; r ¼ r b Þ 2p r b ðD < z < D þ H; r ¼ r b Þ T ¼ Tb _0 ~ TÞ ¼ Q ~ ðD < z < D þ H; r ¼ r sc Þ n  ðksc r 2pr sc

ð9aÞ ð9bÞ

z z ¼ rb ¼

r r ¼ rb 

s ¼

ðT  T gr;1 Þ  kgr Q_ 0



agr  s r 2b

¼ Forb

T Q ðor

r

Q′ 2π rb

n ⋅ k ∇T =

n ⋅ ∇T = 0

n ⋅ ∇T = 0

z (a)

rgr rb D

0

T = Tgr ,∞

r

n ⋅ ∇T = 0 T = Tb

H H gr

n ⋅ ∇T = 0

n ⋅ ∇T = 0

z

(b)

rgr rb D

0

r

n ⋅ ∇T = 0

H H gr

T = Tgr ,∞

n ⋅ k sc ∇ T =

Q′ 2π rsc

n ⋅ ∇T = 0

rsc

ð9cÞ

By introducing suitable dimensionless variables, it is possible to convert the previous equations in corresponding dimensionless ones. With reference to (9a) and (9c) boundary conditions, it is possible to define the following dimensionless variables: 

T = Tgr ,∞

n ⋅ ∇T = 0

H gr

of the thermo-physical properties: ksc ¼ 107  kgr , ðqcÞsc ¼ 105  ðqcÞgr ; in this material, the equation to be solved is again the unsteady Fourier equation without heat generation in cylindrical coordinates. The thermally super-conductive material ensures that at the BHE surface (r = rb) the heat transfer rate is still constant as in the first condition but also the temperature is uniform along the depth, even if it is evolving in time [31].

~ TÞ ¼ ~ n  ðkgr r

0

T sc Þ ð10Þ

n ⋅ ∇T = 0

z (c) Fig. 1. Schematics of the computational domain, with the three different boundary conditions applied at r = rb and at r = rsc: (a) constant heat flux per unit length Q_ 0 ; (b) constant temperature Tb; (c) constant heat flux applied on the superconductive material. Not to scale.

938

A. Priarone, M. Fossa / Applied Thermal Engineering 103 (2016) 934–944 

and the additional quantities r sc ¼ r sc =rb , k ¼ ksc =kgr and ðqcÞ ¼ ðqcÞsc =ðqcÞgr . Worth noticing, the variable r⁄ is the corresponding variable p in the Ingersoll et al. notation for the G and F solutions. The subscript sc indicate the temperature field obtained with the superconductive material boundary condition. With reference to the problem represented by boundary condition (9b), the following dimensionless quantities have been introduced:

z ¼

z rb

r ¼

r rb

agr  s

s ¼

r2b

T T ¼

¼ Forb

ðT  T gr;1 Þ ðT b  T gr;1 Þ

ð11Þ

In the dimensionless 2D axial-section of the geometry, the surface replacing the BHE periphery has radius r⁄ = 1 and length H⁄ = H/rb, while the ground is a cylinder, coaxial with the BHE, with external radius r gr ¼ r gr =rb and length Hgr ¼ Hgr =r b . The following values of the dimensionless BHE length, H⁄ have been considered: 400, 1000, 2000, 3000, 4000. In the study, a fixed value of the D/H ratio equal to 0.04 has been also considered. Defining a dimensionless adiabatic depth D ¼ D=rb , due to the selected H⁄ values this additional parameter resulted to be 16, 40, 80, 120, 160. By means of the dimensionless temperature distribution T Q , T sc or T T obtained through the simulations, one can evaluate by integration the dimensionless mean temperatures T b;Q , T b;sc or T b;T at the boundary surface of the BHE (r = rb) as:

T b;Q ; T b;sc or T b;T ¼

1 H

Z H



T  ðr  ¼ 1Þdz

ð12Þ

Similar integrations, applied to parallel cylindrical surfaces of length H⁄ (and adiabatic depth D⁄) with radius r⁄ > 1 (or equivalently p > 1), allow to evaluate the dimensionless average (along the depth) temperatures T p;Q , T p;sc or T p;T , able to describe the temperature field in the ground domain and also useful for performing spatial superposition of solutions. The values have been calculated for the following values of p: 2, 5, 10, 20, 30, 40, 60, 100, 130, 160, 200, 300, 400, 600, 800, 1000, 2000, 3000. For the imposed temperature boundary condition, one can calculate the average heat transfer rate (per unit length and variable in time) along the BHE depth H⁄:

1 Q_ 0b ¼  H

Z H

3.2. Validation against analytical solution The numerical solution in case of imposed heat flux has been validated in terms of temperature response factor CQ by means of the comparison with the theoretical values described by Eqs. (2)–(4). The validation has been performed for different values of H⁄ and for selected values of the p parameter, namely p = 1, 100, 3000. D/H is set to 0.04. Fig. 2 presents the comparison among numerical and analytical data made for H⁄= 1000, 2000, 4000; data are evaluated only for Forb P 1000, which corresponds, for H⁄= 2000, to about ln(9 FoH) P 6 (the typical bottom end FoH value in g-function representation). The Figure shows that the agreement between numerical and analytical solutions is very good, being the average difference among analytical and numerical solutions equal to 0.05%, 0.25% and 0.95% for p = 1, 100, 1000, respectively; the standard deviation of percentage differences is equal to 0.03%, 0.25% and 0.5% for p = 1, 100, 1000, respectively. 4. Results and discussion Fig. 3 reports a comparison between the temperature transfer functions C evaluated on the BHE (p = 1) pertaining the three different boundary conditions for the case H⁄ = 2000 and for two values of the adiabatic depth, D/H = 0, 0.04. The figure reveals that in the early period, with a threshold Forb which increases with H⁄ (approximately Forb = 1000 for H⁄ = 2000), all the temperature transfer functions C with heat flux boundary condition (with or without superconductive material) overlap. In the same Forb range also the T-boundary condition curves collapse into one, irrespective of the value of the adiabatic depth. Moreover in this early period the temperature transfer functions C at imposed heat flux show lower values compared to the corresponding ones at imposed temperature. The Appendix A reports in tabular form the data presented in Fig. 3 pertaining to no adiabatic depth (D/H = 0) and to the boundary conditions of imposed heat flux and imposed temperature; the Appendix A presents also the same data for H⁄ = 400 and 4000 as a function of the Fourier number and for different values of the dimensionless distance p from the heat source thus providing tables similar to the G and F function ones calculated by Ingersoll et al. for their infinite cylindrical sources.

8

Q_ 0 ðr ¼ 1Þ  dz kgr ðT  T gr;1 Þ

ð13Þ

Finally it is possible to introduce the following temperature response factors according to Eq. (7) definition:

Cb;sc ¼ 2p 

T b;sc

Cb;T

2p ¼ Q_ 0

ð14Þ

b

A period of operation s⁄ = Forb = 10,000,000 have been considered and non-uniform time steps have been adopted in computations [30]. In order to obtain the steady state response factor of the ground, for all the considered H⁄ values and for the three different boundary conditions, a similar set of equations have been solved according to the steady state form of the Fourier equation as a function of the r⁄ and z⁄ variables. Finally, with reference to the boundary condition (9c, superconductive core material), the value of the dimensionless temperature T sc has been numerically calculated for r = rb and its uniformity along z has been verified for different Forb.

6

ΓQ, numerical

Cb;Q ¼ 2p 

T b;Q

D/H = 0.04

7

5

H* = 1000, p = 1 H* = 1000, p = 100 H* = 1000, p = 3000 H* = 2000, p = 1 H* = 2000, p = 100 H* = 2000, p = 3000 H* = 4000, p = 1 H* = 4000, p = 100 H* = 4000, p = 3000

4 3 2 1 0

0

1

2

3

4

5

6

7

8

ΓQ, Claesson Javed Fig. 2. Comparison in terms of temperature transfer functions C between analytical FLS solutions [19] and numerical FCS ones for H⁄ = 1000, 2000, 4000 and p = 1, 100, 3000. Adiabatic depth D/H is set to 0.04.

939

A. Priarone, M. Fossa / Applied Thermal Engineering 103 (2016) 934–944

7 6 5

Γb

4 3

7 6.9 H * = 2000, p = 1 6.8 6.7 6.6 6.5 6.4 6.3 Γb 6.2 6.1 6 5.9 5.8 5.7 5.6 5.5 1.E+05 1.E+06

H * = 2000, p = 1 D/H= 0 Γb,QQ, , D/H =0 T, D/H==00 Γb,T, D/H D/H ==00 Γb,scsc, , D/H Q, = 0.04 Γb,Q, D/HD/H = 0.04 T, D/H= 0.04 = 0.04 Γb,T, D/H Γb,sc, D/H = 0.04 sc, D/H = 0.04

2 1 0 0.001

1

1000

1000000

1.E+08

1.E+09

5.5 5.4 H * = 400, p = 1 5.3 5.2 5.1 5 4.9 Q, ,D/H = 0= 0 Γb 4.8 Γb,Q D/H 4.7 T, D/H = 0= 0 Γ , D/H 4.6 b,T 4.5 sc, ,D/H = 0= 0 Γb,sc D/H 4.4 Q, ,D/H = 0.04 Γb,Q D/H = 0.04 4.3 T, D/H = 0.04 Γ , D/H = 0.04 b,T 4.2 4.1 Γb,sc, D/H = 0.04 sc, D/H= 0.04 4 1.E+05 1.E+06 1.E+07 1.E+08

1.E+09

1E+09

1.E+07

Forb

Forb

(a)

Fig. 3. Comparison between temperature transfer functions obtained with different boundary conditions, with and without adiabatic depth (D/H = 0, 0.04). H⁄ = 2000, p = 1.

Fig. 4a represents a magnification of Fig. 3 at Forb > 105. The presence of a finite adiabatic depth increases the values of each C function and reduces the differences among the results pertaining different boundary conditions. For both the considered D/H values the imposed heat flux condition yields the higher values of the temperature transfer functions (in agreement with [32]). It is also possible to notice that in the late period temperature boundary condition and sc-one provide similar C function trends, with an almost perfect agreement at steady state for D/H = 0.04. Fig. 4b and c depicts again the late and steady state behavior related to dimensionless depths 400 and 4000, respectively. Inspection of Fig. 4 shows that in presence of the adiabatic depth (D/H = 0.04) the temperature transfer functions C related to imposed temperature and to the superconductive material again almost coincide, irrespective of H⁄. On the contrary, in case of no buried depth (D/H = 0), the value of H⁄ plays a role in C function profiles in the late period. In particular it can observed that the differences between CT and Csc increase with the decreasing H⁄. 4.1. Effects of BHE length H⁄ As reported recently also by Cimmino et al. [28], the temperature transfer functions C are usually evaluated at H⁄ref = 2000 which corresponds to a typical borehole length equal to 100 m when BHE diameter is 0.1 m. Eskilson [13] proposed a correction for small variations of H⁄ according to which once a reference solution is available, any other transfer function C can be calculated as:

H Cb;H ðFoH Þ  Cb;Href ðFoH Þ ¼ ln Href

Q, D/H = 0= 0 Γb,Q , D/H T, D/H = Γb,T, D/H0= 0 sc, ,D/H = 0= 0 Γb,sc D/H Q, D/H = 0.04 Γb,Q , D/H = 0.04 T, ,D/H = 0.04 Γb,T D/H = 0.04 = 0.04 Γb,sc, D/H sc, D/H= 0.04

! ð15Þ

The ‘‘true” values of the temperature transfer functions C, numerically obtained by means of Comsol model for lengths H⁄ different from H⁄ref = 2000 have been compared with those calculated using correlation (15). Table 1 reports the values of the maximum relative error introduced by Eq. (15) for the case p = 1 at steady state. The table shows that correlation (15) allows an accurate method for inferring C values when reference ones are available. The error here estimated is in the 1–4% range, depending on the boundary condition. At p = 1 accuracy of correlation (15) demonstrated to be higher at lower H⁄ values. An exception is for short heat sources (H⁄ = 400) with sc-boundary condition: in this case the relative error resulted to be 9%.

Forb

(b) 7.5 7.4 H * = 4000, p = 1 7.3 7.2 7.1 7 6.9 6.8 Γb 6.7 6.6 6.5 6.4 6.3 6.2 6.1 6 1.E+05 1.E+06

Q, ,D/H = 0= 0 Γb,Q D/H T, D/H = Γb,T, D/H 0= 0 sc, ,D/H = 0= 0 D/H Γb,sc Q, ,D/H = 0.04 Γb,Q D/H = 0.04 T, ,D/H = 0.04 Γb,T D/H = 0.04 sc, D/H= 0.04 Γb,sc, D/H = 0.04

1.E+07

1.E+08

1.E+09

Forb

(c) Fig. 4. Comparison between temperature transfer functions obtained with different boundary conditions, with and without adiabatic depth (D/H = 0, 0.04). Forb > 105, p = 1. (a) H⁄ = 2000, (b) H⁄ = 400, (c) H⁄ = 4000.

Eq. (15) is based on a steady state concept of cylindrical shell thermal resistance and hence its validity cannot be extended far from the heat source. The present calculations showed that for p > 1 Eq. (15) cannot be employed without introducing meaningful errors: the relative error increases as Forb decreases and for p increasing values. As an example, the average values of the maximum error induced by Eq. (15) are about 17% and 97% for

940

A. Priarone, M. Fossa / Applied Thermal Engineering 103 (2016) 934–944

Table 1 Comparison between numerically obtained temperature transfer functions C and those calculated using correlation (15). p = 1, steady state values. Boundary condition

Imposed heat flux on the BHE

Imposed temperature on the BHE

Imposed heat flux on the super-conductive material

D/H

0

0.04

0

0.04

0

0.04

p=1 Maximum error induced by Eq. (15) (%) H⁄ with maximum error

1.05 3000

1.29 3000

9.0 400

2.40 3000

3.85 3000

2.73 3000

p = 100 and p = 3000, respectively at H⁄= 400 and again steady state.

4.2. Effect of the boundary condition at BHE surface Fig. 5 shows the effects of the boundary conditions on the temperature transfer functions C as evaluated without adiabatic length (D/H = 0). The figure in particular presents the profile of the ratios Cb,Q/Cb,T, Cb,T/Cb,sc and Cb,Q/Cb,sc as a function of Forb when the dimensionless depth H⁄ = 2000. Fig. 5a refers to the BHE surface (p = 1) whereas Fig. 5b to the case p = 100.

1.8 1.6 1.5

Γb,T/Γb,sc

1.4

Γb ratio

1.3 1.2

Γb,Q/Γb,sc

1.1 1 0.9 0.8

Γb,Q/Γb,T

0.7

4.4. Effect of the adiabatic depth D⁄ and of the distance p from BHE

0.6 0.5 0.001

1

1000

1000000

Forb

(a) 1.8

p = 100, H* = 2000

1.7 1.6 1.5

Γ ratio

1.4 1.3

ΓT/Γsc

1.2 1.1 1 0.9

ΓQ/Γsc

0.8

ΓQ/ΓT

0.7 0.6 0.5 1000

10000

4.3. Effect of the adiabatic depth D⁄ The presence of the buried depth D⁄ affects the temperature transfer functions in the direction of increasing their values, irrespective of the boundary condition. Fig. 6a–c show this behavior in terms of ratios (CD ¼0 =CD ¼0:04 ) as a function of Forb. Data refer to BHE surface (p = 1). As can be noticed from figure inspection, the presence of an adiabatic depth is almost negligible for Forb lower than 10 and more remarkable for short heat sources (H⁄ = 400). The ratios are always very close to unity (apart numerical scatter of data in the first time steps of simulations, i.e. for Forb < 0.1) up to the above threshold dimensionless time: later in time toward infinity the ratios decrease down to 0.96–0.93 when a pure temperature boundary condition is applied and the BHE source is short compared to its radius.

p = 1, H* = 2000

1.7

Concerning Fig. 5a, it can be noticed that in the early period (up to Forb = 103) the Cb,T function is always higher in value than the imposed heat flux ones (up to 1.7 times at heat source activation) that in turn show a similar profile (their ratio is close to the unity). In the late period (Forb > 103) the three Cb functions converge and the related ratios shrink to values close to unity as far as the steady state condition. Fig. 5b, pertaining to p = 100, presents similar trends but the discrepancy between different boundary conditions at the steady condition (and for very high Forb values) increases moving away from the heat source.

100000

1000000

10000000

Forb

(b) Fig. 5. Temperature transfer function ratios as function of dimensionless time: effect of the boundary condition (H⁄ = 2000, D/H = 0). (a) p = 1, (b) p = 100.

Finally the joint effects of adiabatic depth D⁄ and distance p are presented in Fig. 7a–c where the dimensionless height H⁄ is the parameter. Numerical results are here presented as C values at D⁄ = 0 vs the corresponding data at D/H = 0.04 for different distances from the BHE, namely p = 1, 100 and 3000. The Figures refer to super-conductive material, but the displayed trends are similar for the other boundary conditions. As can be noticed Fig. 7a represents in a different form the same results of Fig. 6c; on the other hand Fig. 7b and c show the data for different distances from BHE and clearly suggest that the effect of the adiabatic depth is more noticeable moving toward the steady state period (here represented by the highest values of the C functions) and moving away from the heat source (i.e. at increasing p). This is an important evidence that explains the behavior of multiBHE temperature transfer functions (i.e. the g-function family) in the late period when an adiabatic depth is considered: superposition of temperature fields due to multiple heat sources is driven in this case by either p = 1 or p  1 solutions, thus emphasizing the effect of the D⁄ parameter on the shape of the temperature transfer function. As a preliminary result Fig. 8 reports the temperature response factor pertaining a regular matrix of 64 heat sources having H⁄ = 2000, D⁄ = 0, B⁄ = 100 as obtained by spatial superposition of available solutions. The result of the superposition is the generation of g-functions [13] that in present work notation are referred

941

A. Priarone, M. Fossa / Applied Thermal Engineering 103 (2016) 934–944

8

p= 1

1.02

6 1 0.98

H* = 400 H* = 1000 H* = 2000 H* = 3000 H* = 4000

0.96 0.94 0.92 0.001

1000

H* = 400

2

H* = 2000

H* = 1000 H* = 3000 H* = 4000 0

1

2

3

4

5

6

7

8

Γsc D/H=0

(a)

3

p = 100 2.5

p= 1

Γsc D/H=0.04

D/H = 0)/Γb,T (D/H = 0.04)

3

1000000

1.02

Γb,T

4

0

1

1.04

1 0.98 0.96 0.94 0.92 0.001

H* = 400 H* = 1000 H* = 2000 H* = 3000 H* = 4000

2 1.5

H* = 400

1

H* = 1000 H* = 2000

0.5 0

H* = 3000 H* = 4000 0

0.5

1

1.5

2

1

1000

0.2

p = 3000 0.16

Γsc D/H=0.04

1.02 1

0.12

H* = 400

0.08

H* = 1000

0.98

0.94 0.92 0.001

3

(b)

1000000

p= 1

0.96

2.5

Γsc D/H=0

Forb (b) 1.04

5

1

Forb (a)

Γb,sc (D/H = 0)/Γb,sc (D/H = 0.04)

p=1

7

Γsc D/H=0.04

Γb,Q (D/H = 0)/Γb,Q (D/H = 0.04)

1.04

H* = 2000

0.04

H* = 400 H* = 1000 H* = 2000 H* = 3000 H* = 4000

H* = 3000 H* = 4000

0

0

0.04

0.08

0.12

0.16

0.2

Γsc D/H=0

(c)

1

1000

1000000

Forb (c) Fig. 6. Effect of the adiabatic depth D⁄ on the Cb functions as Forb is varied. (a) Imposed heat flux, (b) Imposed temperature, (c) Super-conductive material.

as Cb,multi TFRs. Interesting to notice in Fig. 8 is the late behavior of the two functions that have been obtained by superposing once CT

Fig. 7. Effect of the adiabatic depth D⁄ and distance p on the Cp functions, H⁄ is the parameter. Super-conductive material (a) p = 1, (b) p = 100, (c) p = 3000.

solutions and CQ ones. The two g-functions have very similar profiles in dimensionless time until ln(9FoH) becomes close to zero: from here on the heat rate based solution shows higher values (around 5%) than the temperature solution as also recently discussed by Cimmino et al. [28].

942

A. Priarone, M. Fossa / Applied Thermal Engineering 103 (2016) 934–944

68

ΓQ superposition

Γmulti or g-function [-]

60

ΓT superposition

52 44 36 28

8x8 matrix of FCS

20

H *=2000, D*=0 p=1, B*=100

12 4

-6

-4

-2

0

2

4

ln(9FoH) Fig. 8. Multiple source temperature transfer functions (g-functions) as obtained from single source solution superposition. Effects of the boundary condition.

5. Conclusions and future work In this paper a numerical model has been developed in order to calculate the temperature response factors pertaining the single BHE when different operating modes in terms of BHE wall boundary conditions are considered. Such conditions are the classic constant heat flux (I) but also the imposed temperature (II) and the imposed heat flux at a superconductive core (III), the latter also satisfying the uniform temperature condition along the heat source depth. The numerical scheme demonstrated to be robust and reliable and able to fit the analytical available solutions (that only

refer to boundary condition I) with great accuracy (average error below 0.35%). New temperature response factors here evaluated include a solution named CT, which refers to the finite length cylindrical source at imposed temperature for which tabulated data have been provided, hence offering an original contribution to the early studies by Ingersoll et al. and Carlslaw and Jaeger. The analysis of the dimensionless temperature profiles in time and space revealed interesting features confirming the possibility to evaluate the temperature response factor at the heat source/ground interface for a single BHE radius to depth ratio and inserting a correction term based on a steady state concept. The results show that the temperature response factor profiles at BHE surface pertaining to conditions (I) and (III) are very close to each other in the early period. Approaching the steady state period condition (I) profiles are always higher (0.80–3.79% at p = 1) than those pertaining to condition (III) which in turn get closer to the temperature condition functions (II) especially when a buried depth is present. Based on above considerations, the imposed heat flux case with uniform temperature (the superconductive core mode) seems to have all the characteristics for spatial and temporal superposition addressed to the calculation of temperature response factors pertaining to multiple BHE arrangements (i.e. g-functions). These aspects will be investigated in a next paper and a careful comparison will be presented with either reference to (I) and (II) condition g-function profiles or the original Eskilson ones. Appendix A. Temperature transfer function values at different dimensionless time and radius for different boundary conditions

CQ (boundary condition: constant heat transfer rate) 400

2000

4000

Forb

p=1

p=2

p=5

p = 10

p = 100

p=1

p=2

p=5

p = 10

p = 100

p=1

p=2

p=5

p = 10

p = 100

0.001 0.005 0.01 0.05 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 1 1.2 1.5 2 2.5 3 4 5 6 8 10 12 15 20 25 30 50 100 500 1000 5000 10,000

0.0337 0.0742 0.1043 0.2242 0.3091 0.4243 0.5035 0.5654 0.6167 0.6611 0.7015 0.7381 0.8025 0.8576 0.9275 1.0218 1.0976 1.1621 1.2671 1.3551 1.4303 1.5524 1.6481 1.7258 1.8222 1.9509 2.0511 2.1331 2.3631 2.6889 3.4301 3.7272 4.3726 4.6003

0 0 0 0.0004 0.0037 0.0196 0.0439 0.0711 0.0988 0.1256 0.1509 0.1750 0.2200 0.2612 0.3165 0.3954 0.4617 0.5193 0.6151 0.6963 0.7660 0.8812 0.9730 1.0484 1.1423 1.2683 1.3669 1.4478 1.6757 1.9996 2.7393 3.0362 3.6815 3.9092

0 0 0 0 0 0 0 0 0 0 0.0001 0.0003 0.0009 0.0019 0.0046 0.0120 0.0225 0.0352 0.0649 0.0956 0.1265 0.1860 0.2415 0.2925 0.3609 0.4582 0.5392 0.6082 0.8121 1.1143 1.8368 2.1320 2.7755 3.0031

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0003 0.0010 0.0024 0.0074 0.0152 0.0254 0.0437 0.0784 0.1152 0.1517 0.2829 0.5153 1.1762 1.4645 2.1013 2.3283

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0005 0.0108 0.2181 0.3750

0.0343 0.0749 0.1042 0.2239 0.3090 0.4246 0.5040 0.5661 0.6175 0.6619 0.7024 0.7391 0.8036 0.8590 0.9292 1.0241 1.1004 1.1651 1.2711 1.3590 1.4349 1.5581 1.6551 1.7338 1.8313 1.9617 2.0637 2.1474 2.3834 2.7199 3.5013 3.8372 4.6126 4.9267

0 0 0 0.0004 0.0037 0.0197 0.0440 0.0712 0.0989 0.1259 0.1513 0.1754 0.2206 0.2618 0.3174 0.3967 0.4634 0.5212 0.6179 0.6992 0.7696 0.8857 0.9787 1.0550 1.1501 1.2778 1.3781 1.4607 1.6945 2.0291 2.8090 3.1446 3.9199 4.2340

0 0 0 0 0 0 0 0 0 0 0.0001 0.0003 0.0009 0.0019 0.0046 0.0120 0.0226 0.0354 0.0652 0.0963 0.1274 0.1875 0.2434 0.2950 0.3643 0.4630 0.5453 0.6157 0.8247 1.1368 1.8992 2.2325 3.0060 3.3199

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0003 0.0011 0.0024 0.0075 0.0154 0.0256 0.0442 0.0794 0.1167 0.1539 0.2880 0.5277 1.2267 1.5515 2.3180 2.6312

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0006 0.0123 0.2675 0.4928

0.0349 0.0743 0.1031 0.2233 0.3067 0.4212 0.5031 0.5669 0.6190 0.6640 0.7036 0.7389 0.8022 0.8574 0.9288 1.0255 1.1030 1.1685 1.2745 1.3607 1.4317 1.5534 1.6514 1.7341 1.8346 1.9631 2.0654 2.1505 2.3890 2.7142 3.5126 3.8472 4.6400 4.9722

0 0 0 0.0005 0.0041 0.0205 0.0441 0.0707 0.0982 0.1250 0.1508 0.1754 0.2208 0.2618 0.3175 0.3974 0.4647 0.5230 0.6200 0.7000 0.7671 0.8822 0.9757 1.0553 1.1529 1.2790 1.3797 1.4636 1.6998 2.0234 2.8201 3.1545 3.9471 4.2793

0 0 0 0 0 0 0 0 0 0.0001 0.0001 0.0003 0.0009 0.0019 0.0046 0.0120 0.0224 0.0351 0.0646 0.0963 0.1286 0.1887 0.2442 0.2955 0.3649 0.4634 0.5462 0.6174 0.8281 1.1324 1.9093 2.2415 3.0322 3.3642

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0001 0.0003 0.0010 0.0024 0.0075 0.0155 0.0256 0.0440 0.0796 0.1169 0.1540 0.2883 0.5277 1.2350 1.5592 2.3426 2.6737

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0007 0.0122 0.2739 0.5087

943

A. Priarone, M. Fossa / Applied Thermal Engineering 103 (2016) 934–944 (continued)

CQ (boundary condition: constant heat transfer rate) 400

2000

4000

Forb

p=1

p=2

p=5

p = 10

p = 100

p=1

p=2

p=5

p = 10

p = 100

p=1

p=2

p=5

p = 10

p = 100

25,000 100,000 250,000 1,000,000 2,000,000 5,000,000 1.00E+07 Steady state

4.8166 4.9584 4.9855 4.9971 4.9983 4.9991 4.9993 4.9993

4.1255 4.2673 4.2944 4.3060 4.3072 4.3080 4.3082 4.3082

3.2193 3.3610 3.3881 3.3997 3.4010 3.4017 3.4020 3.4020

2.5442 2.6859 2.7130 2.7247 2.7259 2.7266 2.7269 2.7269

0.5598 0.6940 0.7206 0.7322 0.7334 0.7341 0.7344 0.7344

5.3316 5.9000 6.2022 6.4925 6.5536 6.5855 6.5999 6.6056

4.6388 5.2072 5.5095 5.7997 5.8608 5.8928 5.9072 5.9129

3.7246 4.2930 4.5952 4.8855 4.9466 4.9785 4.9929 4.9986

3.0354 3.6035 3.9057 4.1960 4.2571 4.2890 4.3034 4.3091

0.8372 1.3736 1.6712 1.9598 2.0208 2.0528 2.0671 2.0728

5.4000 6.0284 6.4050 6.8948 7.0652 7.2102 7.2649 7.3023

4.7071 5.3354 5.7121 6.2019 6.3723 6.5173 6.5720 6.6093

3.7919 4.4202 4.7968 5.2866 5.4570 5.6020 5.6567 5.6940

3.1009 3.7289 4.1055 4.5953 4.7657 4.9107 4.9654 5.0027

0.8753 1.4691 1.8396 2.3265 2.4967 2.6415 2.6962 2.7335

CQ (boundary condition: constant heat transfer rate) Forb

400 p=1

p=2

p=5

p = 10

p = 100

2000 p=1

p=2

p=5

p = 10

p = 100

4000 p=1

p=2

p=5

p = 10

p = 100

0.001 0.005 0.01 0.05 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 1 1.2 1.5 2 2.5 3 4 5 6 8 10 12 15 20 25 30 50 100 500 1000 5000 10,000 25,000 100,000 250,000 1,000,000 2,000,000 5,000,000 1.00E+07 Steady state

0.0562 0.1188 0.1616 0.3336 0.4410 0.5818 0.6724 0.7402 0.8010 0.8554 0.9016 0.9436 1.0138 1.0708 1.1415 1.2326 1.3128 1.3819 1.4882 1.5719 1.6423 1.7528 1.8383 1.9067 1.9932 2.1030 2.2014 2.2813 2.4993 2.7728 3.4196 3.6818 4.1963 4.3698 4.5161 4.6183 4.6427 4.6521 4.6531 4.6537 4.6539 4.6652

0 0 0 0.0004 0.0079 0.0474 0.0954 0.1412 0.1845 0.2247 0.2615 0.2961 0.3573 0.4092 0.4752 0.5624 0.6391 0.7057 0.8097 0.8922 0.9620 1.0719 1.1572 1.2256 1.3120 1.4221 1.5206 1.6006 1.8194 2.0940 2.7438 3.0073 3.5242 3.6985 3.8455 3.9482 3.9727 3.9821 3.9832 3.9837 3.9840 3.9937

0 0 0 0 0 0 0 0 0 0.0001 0.0003 0.0007 0.0022 0.0048 0.0109 0.0262 0.0454 0.0669 0.1126 0.1577 0.2006 0.2779 0.3443 0.4015 0.4756 0.5754 0.6630 0.7364 0.9454 1.2171 1.8686 2.1344 2.6565 2.8326 2.9811 3.0847 3.1095 3.1190 3.1201 3.1206 3.1208 3.1284

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0002 0.0008 0.0024 0.0052 0.0145 0.0281 0.0445 0.0722 0.1206 0.1674 0.2119 0.3630 0.6005 1.2258 1.4912 2.0173 2.1951 2.3451 2.4497 2.4746 2.4843 2.4854 2.4859 2.4861 2.4922

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0008 0.0136 0.2263 0.3704 0.5149 0.6213 0.6471 0.6571 0.6582 0.6588 0.6590 0.6606

0.0567 0.1197 0.1658 0.3358 0.4444 0.5819 0.6751 0.7485 0.8073 0.8588 0.9044 0.9450 1.0136 1.0711 1.1417 1.2469 1.3285 1.3950 1.4985 1.5851 1.6581 1.7736 1.8632 1.9360 2.0348 2.1613 2.2569 2.3331 2.5528 2.8837 3.6050 3.9275 4.6495 4.9377 5.3300 5.8335 6.0913 6.3155 6.3664 6.4017 6.4148 6.4195

0 0 0 0.0004 0.0081 0.0478 0.0960 0.1427 0.1857 0.2254 0.2620 0.2958 0.3556 0.4075 0.4732 0.5709 0.6486 0.7128 0.8141 0.8990 0.9708 1.0851 1.1741 1.2465 1.3448 1.4709 1.5663 1.6424 1.8620 2.1929 2.9145 3.2372 3.9597 4.2481 4.6407 5.1445 5.4025 5.6268 5.6777 5.7130 5.7261 5.7308

0 0 0 0 0 0 0 0 0 0.0001 0.0003 0.0007 0.0022 0.0048 0.0110 0.0265 0.0463 0.0684 0.1147 0.1603 0.2035 0.2818 0.3499 0.4093 0.4900 0.5993 0.6863 0.7578 0.9675 1.2883 2.0064 2.3291 3.0525 3.3413 3.7346 4.2392 4.4977 4.7224 4.7734 4.8088 4.8219 4.8266

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0001 0.0007 0.0024 0.0053 0.0151 0.0292 0.0461 0.0742 0.1240 0.1731 0.2194 0.3750 0.6411 1.3291 1.6482 2.3703 2.6596 3.0534 3.5591 3.8180 4.0431 4.0942 4.1296 4.1428 4.1475

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0009 0.0161 0.3019 0.5301 0.8789 1.3754 1.6362 1.8638 1.9155 1.9514 1.9646 1.9694

0.0568 0.1228 0.1673 0.3385 0.4446 0.5847 0.6764 0.7487 0.8083 0.8601 0.9052 0.9445 1.0132 1.0725 1.1464 1.2448 1.3259 1.3935 1.5011 1.5835 1.6580 1.7798 1.8721 1.9457 2.0376 2.1677 2.2672 2.3468 2.5757 2.8828 3.6374 3.9567 4.7070 5.0317 5.4440 6.0302 6.3844 6.8126 6.9587 7.0743 7.1201 7.1515

0 0 0 0.0005 0.0081 0.0479 0.0962 0.1429 0.1861 0.2259 0.2623 0.2956 0.3554 0.4083 0.4761 0.5687 0.6460 0.7111 0.8160 0.8971 0.9704 1.0906 1.1821 1.2553 1.3469 1.4765 1.5757 1.6552 1.8839 2.1908 2.9455 3.2650 4.0155 4.3403 4.7528 5.3391 5.6934 6.1218 6.2680 6.3836 6.4293 6.4607

0 0 0 0 0 0 0 0 0 0.0001 0.0003 0.0007 0.0022 0.0049 0.0112 0.0268 0.0467 0.0687 0.1150 0.1606 0.2043 0.2837 0.3528 0.4127 0.4905 0.6027 0.6924 0.7663 0.9842 1.2843 2.0338 2.3530 3.1037 3.4288 3.8416 4.4285 4.7831 5.2118 5.3581 5.4738 5.5196 5.5509

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0001 0.0007 0.0023 0.0052 0.0150 0.0291 0.0461 0.0746 0.1244 0.1739 0.2209 0.3800 0.6384 1.3517 1.6677 2.4164 2.7415 3.1546 3.7420 4.0970 4.5260 4.6725 4.7883 4.8341 4.8654

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0010 0.0165 0.3133 0.5597 0.9316 1.5040 1.8585 2.2899 2.4373 2.5539 2.6000 2.6311

References [1] B.M.S. Giambastiani, F. Tinti, D. Mendrinos, M. Mastrocicco, Energy performance strategies for the large scale introduction of geothermal energy in residential and industrial buildings: the GEO.POWER project, Energy Policy 65 (2014) 315–322. [2] A. García, E. Santoyo, Paredes G. Espinosa, I. Hernández, H. Gutiérrez, Estimation of temperatures in geothermal wells during circulation and shutin in the presence of lost circulation, Transp. Porous Media 33 (1998) 103–127. [3] N. Molina-Giraldo, P. Blum, K. Zhu, P. Bayer, Z. Fang, A moving finite line source model to simulate borehole heat exchangers with groundwater advection, Int. J. Therm. Sci. 50 (2011) 2506–2513. [4] H.S. Carslaw, J.C. Jaeger, Conduction of Heat in Solids, Claremore Press, Oxford, U.K., 1947. [5] L.R. Ingersoll, O.J. Zobel, A.C. Ingersoll, Heat Conduction with Engineering, Geological, and other Applications, McGraw-Hill, New York, 1954.

[6] M. Philippe, M. Bernier, D. Marchio, Validity ranges of three analytical solutions to heat transfer in the vicinity of single boreholes, Geothermics 38 (2009) 407–413. [7] L. Lamarche, Short-term behavior of classical analytic solutions for the design of ground-source heat pumps, Renewable Energy 57 (2013) 171–180. [8] P. Mogensen, Fluid to duct wall heat transfer in duct system heat storage, in: International Conference on Surface Heat Storage in theory and Practice, Stockholm, Sweden, 1983, pp. 652–657. [9] S. Gehlin, Thermal Response Test, in-situ Measurements of Thermal Properties in Hard Rock Licentiate thesis, Lulea University, 1998. [10] D. Pahud, B. Matthey, Comparison of the thermal performance of double Upipe borehole heat exchangers measured in situ, Energy Build. 33 (2001) 503– 507. [11] J. Acuna, P. Mogensen, B. Palm, Distributed thermal response test on a U-pipe borehole heat exchanger, in: The 11th Int. Conf. on Energy Storage EFFSTOCK, Stockholm, 2009.

944

A. Priarone, M. Fossa / Applied Thermal Engineering 103 (2016) 934–944

[12] S.P. Kavanaugh, K. Rafferty, Ground-Source Heat Pumps: Design of Geothermal Systems for Commercial and Institutional Buildings, ASHRAE, Atlanta, Georgia, 1997. [13] P. Eskilson, Thermal Analysis of Heat Extraction Boreholes Ph.D. Thesis, Lund University of Technology, Sweden, 1987. [14] L. Lamarche, B. Beauchamp, A new contribution to the finite line-source model for geothermal boreholes, Energy Build. 39 (2007) 188–198. [15] H.Y. Zeng, N.R. Diao, Z.H. Fang, A finite line-source model for boreholes in geothermal heat exchangers, Heat Transfer-Asian Res. 31 (2002) 558–567. [16] T.V. Bandos, A. Montero, E. Fernández, J.L.G. Santander, J.M. Isidro, J. Pérez, P.J. Fernández de Córdoba, J.F. Urchueguí, Finite line-source model for borehole heat exchangers: effect of vertical temperature variations, Geothermics 38 (2009) 263–270. [17] M. Fossa, The temperature penalty approach to the design of borehole heat exchangers for heat pump applications, Energy Build. 43 (2011) 1473–1479. [18] M. Fossa, D. Rolando, Fully analytical FLS solution for fast calculation of temperature response factors in geothermal heat pump borefield design, in: 11th International Energy Agency Heat Pump Conference, Montreal, May 12– 16, 2014, 2014. [19] J. Claesson, S. Javed, An analytical method to calculate borehole fluid temperatures for time-scales from minutes to decades, ASHRAE Trans. 117 (2) (2011) 279–288. [20] G. Hellstrom, Duct Ground Heat Storage Model, Manual for Computer Code, Department of Mathematical Physics, Lund Institute of Technology, Lund, Sweden, 1989. [21] G. Hellström, B. Sanner, Earth Energy Designer: Software for Dimensioning of Deep Boreholes for Heat Extraction, Department of Mathematical Physics, Lund University, Sweden, 1994.

[22] J.D. Spitler, GLHEPRO – a design tool for commercial building ground loop heat exchangers, in: Proceedings of the Fourth International Heat Pump in Cold Climates Conference, Aylmer, Québec, 2000. [23] E. Zanchini, S. Lazzari, New g-functions for the hourly simulation of double Utube borehole heat exchanger fields, Energy 70 (2014) 444–455. [24] C. Yavuzturk, J.D. Spitler, A short time step response factor model for vertical ground loop heat exchangers, ASHRAE Trans. 105 (1999) 475–485. [25] M. Bernier, P. Pinel, R. Labib, R. Paillot, A multiple load aggregation algorithm for annual hourly simulations of GCHP systems, Int. J. Heating, Ventilating, AirConditioning and Refrigeration Research 10 (2004) 471–488. [26] F. Sheriff, Génération de facteurs de réponse pour champs de puits géothermiques verticaux M.Sc.A., Thesis, Département de Génie Mécanique. École Polytechnique de Montréal, 2007. [27] O. Cauret, M. Bernier, Experimental validation of an underground compact collector model, in: Proc. 11th Int. Conf. Effstock Stockholm, 2009. [28] M. Cimmino, M. Bernier, F. Adams, A contribution towards the determination of g-functions using the finite line source, Appl. Therm. Eng. 51 (2013) 1–2. [29] COMSOL AB, COMSOL Multiphysics Modeling Guide Version 3.5a, COMSOL AB, Burlington MA, 2008. [30] A. Priarone, M. Fossa, Modelling the ground volume for numerically generating single borehole heat exchanger response factors according to the cylindrical source approach, Geothermics 58 (2015) 32–38. [31] P. Monzo, J. Acuna, A novel numerical model for the thermal response of borehole heat exchanger fields, in: 11th IEA Heat Pump Conference, Montréal (Québec), Canada, 2014. [32] P. Monzo, P. Mogensen, J. Acuna, F. Ruiz-Calvo, C. Montagud, A novel numerical approach for imposing a temperature boundary condition at the borehole wall in borehole fields, Geothermics 56 (2015) 35–44.