A new numerical method to determine isothermal g-functions of borehole heat exchanger fields

Geothermics 77 (2019) 278–287

Contents lists available at ScienceDirect

Geothermics journal homepage: www.elsevier.com/locate/geothermics

A new numerical method to determine isothermal g-functions of borehole heat exchanger fields

T

Claudia Naldi , Enzo Zanchini ⁎

Department of Industrial Engineering (DIN), Alma Mater Studiorum University of Bologna, Viale Risorgimento 2, 40136 Bologna, Italy

ARTICLE INFO

ABSTRACT

Keywords: Ground-coupled heat pumps Borehole heat exchangers Thermal response factor G-function Isothermal boundary condition Numerical method

An iterative numerical method to determine g-functions of borehole fields with the boundary condition of uniform temperature and time-constant surface-averaged heat flux is proposed. The method employs boundary conditions of uniform time-dependent temperature that converge to time-constant surface-averaged heat flux. Two simulations are sufficient for technical purposes, while additional simulations yield very precise g-functions with the desired boundary condition. The method is applied to analyze the overestimation of the g-function yielded by the condition of uniform and constant heat flux and the effects of the buried depth and of the spacing between boreholes on the g-function, for a 3 × 2 borehole field.

1. Introduction Ground-source heat pumps (GSHPs) are becoming an important technology for building heating and cooling. According to Self et al. (2013), GSHPs are not only the most efficient, but also the most convenient system for building heating in Canada. Lund and Boyd, (2016) stated that the annual use of thermal energy produced by GSHPs increased from 24.31 to 90.79 TW h in the decade 2005-2015. Groundcoupled heat pumps (GCHPs) are the most applied GSHP systems. They employ vertical ground heat exchangers, often called borehole heat exchangers (BHEs), having a usual length between 50 and 200 m and a usual diameter of about 15 cm. The simplest method for the design of a BHE field is that developed by Kavanaugh and Rafferty (1997) and recommended in (ASHRAE, 2007 Handbook). A more accurate design and the dynamic simulation of a BHE field can be performed by employing dimensionless thermal response factors, called g-functions, that are often embedded in simulation codes, such as EED (Hellström and Sanner, 1994), GHLEPRO (Spitler, 2000), EnergyPlus (Fisher et al., 2006), eQuest (Liu and Hellström, 2006). G-functions are employed also in research papers on the simulation of GCHP systems, such as (Naldi and Zanchini, 2017, 2018). A g-function is a dimensionless thermal response factor that yields the temperature change at the surface between boreholes and ground produced by a time-constant total heat flux applied to the ground at that surface. With a suitable definition of dimensionless temperature, the g-function of a BHE field is the dimensionless temperature, averaged

⁎

along the BHE length, produced by a time-constant total heat flux applied by the surface of the BHE field to the ground. For a single BHE, gfunctions are often determined not only on the BHE surface, but also on cylindrical surfaces placed at different radial distances from the BHE axis. This is because, under the boundary condition of uniform and constant heat flux at the surface of the BHEs, the g-function for a BHE field can be calculated by the superposition of the effects of the single boreholes. Namely, the temperature change at the surface of a BHE can be determined as the sum of the changes produced, on that surface, by the BHE itself and by all the other BHEs of the field, each considered at its radial distance from the BHE under exam. The simplest method to obtain the g-function for a BHE field is to apply the Finite Line Source (FLS) scheme and the superposition of the effects of the single boreholes. In the FLS scheme, each BHE is considered as a line source with finite length that delivers to the surrounding semi-infinite solid medium a uniform and constant heat flux. The analytical solution of the temperature field produced by a FLS was presented, in form of an integral, by Claesson and Eskilson (1987). The solution was deduced also by Zeng et al. (2002), who pointed out that the use of the solution at the middle of the length of the BHE yields an overestimation of the temperature field. They suggested to employ the expression of the mean value of the temperature along the BHE length, in the form of a double integral. Simpler forms of the FLS solution averaged along the BHE length, in the form of a single integral, were developed by Lamarche and Beauchamp (2007) and by Bandos et al. (2009). Extended solutions, which hold also for BHEs whose top is placed at a buried depth Bd under the ground surface, were presented

Corresponding author. E-mail address: [email protected] (C. Naldi).

https://doi.org/10.1016/j.geothermics.2018.10.007 Received 10 May 2018; Received in revised form 26 September 2018; Accepted 7 October 2018 0375-6505/ © 2018 Elsevier Ltd. All rights reserved.

Geothermics 77 (2019) 278–287

C. Naldi, E. Zanchini

Nomenclature

z

BC Bd BHE D FCS FLS g GCHP GSHP H HCM k p q ql r S T Tg x y

Greek symbols

Boundary condition Borehole buried depth (m) Borehole Heat Exchanger Borehole diameter (m) Finite Cylindrical Source Finite Line Source g-function Ground-Coupled Heat Pump Ground-Source Heat Pump Borehole length (m) High-Conductive-Material method Ground thermal conductivity (W m−1 K−1) Order of accuracy of the discretization scheme Heat flux per unit area (W/m2) Heat flux per unit length (W/m) Grid refinement factor Borehole spacing (m) Temperature (K) Undisturbed ground temperature (K) First coordinate axis Second coordinate axis

α τ

Third coordinate axis

Ground thermal diffusivity (m2/s) Time (s)

Superscripts * –

Dimensionless Averaged on the borehole surface

Subscripts a b CB ees h q rh T

Central boreholes Lateral boreholes Determined by Cimmino and Bernier Estimate of the exact solution Obtained with the fine mesh Isoflux Obtained with the coarse mesh Isothermal

technical applications, where BHEs are feed in parallel with the same inlet temperature, the original Eskilson’s boundary condition has been considered in recent researches. By an analytical method based on the FLS solution, Cimmino et al. (2013) determined g-functions with a boundary condition intermediate between uniform and constant wall heat flux and Eskilson’s condition. They imposed a uniform heat flux along each BHE and equal temperatures for all BHEs. They found results with small differences with respect to the usual boundary condition of uniform and constant wall heat flux. Cimmino and Bernier (2014) developed an improved semianalytical method to determine g-functions with Eskilson’s boundary condition: BHEs are divided into segments releasing different time-dependent heat fluxes per unit length, such that the surface temperature of the BHE field is uniform at each instant of time. The analytical FLS solution is adopted for each segment. The authors compared the results obtained by this method with those obtained by applying the finitedifference numerical method recommended by Eskilson (1987), and found a good agreement between these methods. An improvement of the method employed by Cimmino and Bernier (2014), that takes into account the effects of the BHE thermal resistance on the g-functions, was developed by Cimmino (2015). A technique to reduce the computation time in the application of the method of Cimmino and Bernier (2014), by reducing the number of numerical computations of the FLS solutions, was presented by Cimmino (2018). Lamarche (2017) proposed a modification of the method of Cimmino and Bernier (2014) that allows reducing the number of segments for each borehole. Numerical methods to determine g-functions with Eskilson’s boundary condition were presented by Priarone and Fossa (2016) for a single BHE, and by Monzó et al. (2015) for BHE fields. In (Priarone and Fossa, 2016) the uniform temperature along the BHE was obtained by imposing a uniform and constant heat flux at the inner surface of a superconductive layer having as external radius the BHE radius. In (Monzó et al., 2015) the BHEs were supposed to be filled with a superconductive material, which also interconnects the BHEs and ensures a uniform temperature of the external surface of the BHE field. The authors found a fair agreement between their g-function for a 3 × 2 BHE field and that determined by Cimmino and Bernier (2014). Puttige et al. (2016) and Monzó et al. (2018) proposed an improvement of the numerical model employed by Monzó et al. (2015). The new model takes into account the BHE thermal resistance by introducing a thin resistive

by Costes and Peysson (2008) and by Claesson and Javed (2011). A simple approximate form of the FLS solution, based on regression analysis of many solutions, was proposed by Fossa (2011). Another scheme, suitable to determine g-functions by numerical methods, is to consider each BHE as a cylinder with diameter D and length H that delivers a uniform and constant power per unit length to the surrounding solid medium. Zanchini and Lazzari (2013, 2014) employed this scheme, that they called Finite Cylindrical-Source, FCS, to determine polynomial expressions of the g-functions for a single BHE, for several values of the ratio H/D and of the dimensionless distance r/D from the BHE axis. Zanchini and Lazzari (2014) presented also polynomial expressions of new g-functions, where the BHE grout is included in computations. Monzó et al., 2013 calculated numerically, by the FCS scheme, the g-functions for a square field of 64 BHEs with three different values of the buried depth, and compared the results obtained by these simulations with those determined by the FLS solution and spatial superposition of effects and with those yielded by the g-functions embedded in EED. They found a good agreement between the numerical results and the FLS solution, and higher values of the numerical gfunctions with respect to those obtained through EED. Numerical calculations of the g-functions with a different FCS scheme were performed by Eskilson (1987). The boundary condition adopted by Eskilson is a uniform temperature and a time-constant total heat flux at the surface between the BHE field and the ground. The gfunctions with this boundary condition were evaluated for many fields of long BHEs, with H/D = 1000. The results, reported in graphical form, were obtained by finite-difference numerical simulations. Lamarche and Beauchamp (2007) and Fossa et al. (2009) compared gfunctions obtained by the FLS model and the superposition of effects with the g-functions obtained by Eskilson (1987) for some BHE fields. They found that the g-functions obtained by the FLS solution slightly overestimate Eskilson’s g-functions. Applying the boundary condition employed by Eskilson to determine g-functions is much more difficult than applying the condition of uniform and constant wall heat flux. In fact, the time dependent heat flux distribution that yields a uniform surface temperature and a timeconstant total heat flux is unknown. Moreover, the superposition of the effects of single BHEs does not hold with Eskilson’s boundary condition, so that the g-functions must be calculated for every BHE field. However, since Eskilson’s boundary condition seems closer to that occurring in 279

Geothermics 77 (2019) 278–287

C. Naldi, E. Zanchini

layer between the superconductive material and the ground. In this paper, we propose a new numerical method to obtain the gfunction of a BHE field with Eskilson’s boundary condition, that we will call isothermal g-function. In the new method, the uniformity of the surface temperature of the borehole field is imposed as a boundary condition, and a time-constant total heat flux at the surface of the borehole field is obtained by a recursive procedure, as follows. A simulation of the BHE field with the boundary condition of dimensionless wall heat flux per unit length equal to 1 is performed and the g-function with uniform heat flux is obtained. This g-function is imposed as a boundary condition at the surface of the BHE field, and the time evolution of the dimensionless heat flux per unit length, averaged along the whole BHE field, is determined. The inverse of the time-averaged value of this heat flux, from the initial instant to the current time instant, is used as a correction factor to determine the first-trial isothermal gfunction. This is used as a boundary condition to determine a new dimensionless heat flux and, by applying the new correction factor, the second-trial isothermal g-function. The method can be easily implemented, because the computational domain includes only the ground, with cylindrical holes that represent the BHEs. The surface of the BHE field is perfectly isothermal, because a uniform time-dependent temperature of this surface is imposed as a boundary condition. Therefore, a high accuracy is obtained when a convergence to a timeconstant length-averaged linear heat flux from the BHE field to the ground is reached. After few iterations, the maximum difference between the length-averaged linear heat flux and the desired constant value, and that between the isothermal g-functions obtained in successive iterations become extremely small, so that the method is very precise. The first-trial isothermal g-function, obtained by two simulations, is already sufficiently precise for technical purposes. The method can be easily modified to take into account the BHE thermal resistance, by introducing a thin resistive layer at the surface of the BHEs.

Fig. 2. Computational domain and selected mesh, field with buried depth 4 m and borehole spacing 7.5 m.

Three 3 × 2 BHE fields, with D = 15 cm and H = 150 m are considered: the first with Bd = 4 m and S = 7.5 m, the second with Bd = 1.5 m and S = 7.5 m, and the third with Bd = 1.5 m and S = 6 m. Thanks to the double symmetry of the 3 × 2 BHE field, only a quarter of the field is simulated. A quarter of cylinder, comprising one borehole and a half, is adopted as 3D computational domain. The geometry of the domain and the mesh selected are shown in Fig. 2, which refers to the BHE field with Bd = 4 m and S = 7.5 m. The computational domain is chosen big enough to yield results independent of the domain size, as will be demonstrated in Section 2.1. The outer cylindrical surface and the bottom surface of the domain are considered as adiabatic, while a boundary condition (BC) of symmetry is used for the vertical symmetry planes. The upper surface, namely the horizontal ground surface, is set as isothermal, with a constant temperature equal to the undisturbed ground temperature. The effects of variations over time of the external air temperature are neglected, as is usual for long boreholes. The top and the bottom of the BHEs are considered as adiabatic. The vertical surface of the boreholes is set as isothermal, with time-dependent temperature values obtained by means of a recursive procedure, as will be soon explained. The initial condition is that of uniform temperature equal to the undisturbed ground temperature, for the whole domain. The recursive procedure employed to evaluate the isothermal gfunctions requires a preliminary simulation with the BC of a uniform

2. Numerical method The method presented in this Section can be used to evaluate numerically the isothermal g-function for a BHE field of any geometry. In this paper, three isothermal g-functions are computed with reference to a 3 × 2 borehole field, i.e. a field with 3 in-line BHEs placed in 2 rows. Each borehole is considered as a finite cylindrical source with length H and diameter D, buried at a depth Bd from the ground surface and placed at a distance S from the adjoining BHEs of the field. The ground is considered as a semi-infinite solid medium, with constant thermophysical properties, without groundwater movement. A horizontal view and a vertical section of a 3 × 2 BHE field are reported in Fig. 1, to illustrate the geometrical parameters defined above.

Fig. 1. Horizontal view (left) and vertical section of one row (right) of a 3 × 2 borehole field. 280

Geothermics 77 (2019) 278–287

C. Naldi, E. Zanchini

and time-constant heat flux per unit area, q, at the BHE-ground surface:

q=

k ( T n)|BHE

ground

q = l , D

discretization. The procedure employed to select the mesh and to check the domain-extension independence of the simulation results is reported in Section 2.1.

(1)

where k is the ground thermal conductivity, T is temperature, n is the outward unit normal and ql is the heat flux per unit length applied to the BHE-ground boundary. The differential equation of transient heat conduction to be solved is

T

2 T,

=

2.1. Selected mesh and domain-extension independence The selections of the mesh and of the computational-domain extension are performed with reference to the isoflux g-function of the 3 × 2 field with buried depth 4 m and BHE spacing 7.5 m, with dimensionless geometric parameters H* = H/D = 1000, S* = S/D = 50, B*d = Bd/D = 26.667. The isoflux g-function of the field is calculated by employing two uniformly-refined meshes in order to evaluate the estimate of the exact solution through Richardson extrapolation, according to Richardson (1911, 1927) and Oberkampf and Roy (2010). An unstructured mesh with 289 203 tetrahedral elements is chosen, with minimum dimensionless element size 0.79 in correspondence of the BHEs surface, and maximum dimensionless element size 500 at the external boundaries of the domain. The resulting g-function is denoted by grh. A finer unstructured mesh with a regular refinement factor r = 2 is then built in COMSOL by splitting every tetrahedral element of the previous mesh in 8 smaller elements. The splitting is obtained by connecting the midpoints of each element edge. The finer mesh has 2 313 624 elements. The g-function determined with this mesh is denoted by gh. The extrapolated estimate of the exact solution, gees, is

(2)

where τ is time and α is the ground thermal diffusivity. By introducing the following dimensionless quantities *=D ,

x* =

x , D

y* =

y , D

z* =

z , D

*=

D2

,

T* = k

T

Tg ql

,

(3)

where Tg is the undisturbed ground temperature, one can rewrite Eq. (2) in the dimensionless form

T* = *

*2 T *.

(4)

The dimensionless BC at the BHE-ground interface is

q* =

( *T * n)|BHE

ground

=

1

.

(5)

Since the dimensionless length of the BHE circumference is π, the dimensionless heat flux per unit length is

ql* = q* = 1.

gees = gh +

(6)

The dimensionless BC on the horizontal surface of the ground and the dimensionless initial condition are, respectively,

T * (x *, y *, 0, *) = 0,

T * (x *, y *, z *, 0) = 0.

gh

grh

rp

1

,

(8)

where p is the order of accuracy of the discretization scheme, equal to 3 because the polynomial order of the tetrahedral mesh is two. A comparison between grh, gh and gees is shown in Fig. 3, where the small differences are highlighted by the zoom on the graph. Values of the g-functions grh, gh and gees for selected values of log10 (τ*) are reported in Table 1, together with the relative discrepancy between grh and gees, and that between gh and gees. Computations are performed by a PC with Intel Core i7-6700 K 4.0 GHz, RAM 64 GB. The simulation time is about 14 min with the coarse mesh and about 10 h and 30 min with the fine one. Due to the long simulation time required by the fine mesh, alternative meshes yielding precise results with a shorter computation time are sought. Simulations are performed with two new unstructured meshes: Mesh 1 with 965 056 tetrahedral elements, minimum dimensionless element size 0.64 and maximum one 500; Mesh 2 with 1 350 848 tetrahedral elements, minimum dimensionless element size 0.34 and maximum one 250. The latter is obtained by introducing a local refinement in the neighborhood of the boreholes. The computation time is about 2 h and 10 min with Mesh 1 and about 1 h and 30 min with Mesh 2. The gfunctions obtained by these meshes, g1 and g2, are compared with the

(7)

The dimensionless Eqs. (4,5,7) are solved by means of the software COMSOL Multiphysics, which employs the finite element method. The computational domain (quarter of cylinder) has dimensionless radius 2500 and dimensionless height 2500 + B*d, with dimensionless buried depth B*d = Bd /D. An unstructured mesh with tetrahedral elements is adopted, with finer elements in proximity of the boreholes and coarser elements towards the external boundaries. The selected mesh, which is shown in Fig. 2 for the field with Bd = 4 m and S = 7.5 m, consists of about 1 350 000 elements. The simulation time interval is 10−6 ≤ τ* ≤ 106, namely -6 ≤ log10 (τ*) ≤ 6 in logarithmic scale, divided in time steps of log10 (τ*) = 0.05. With a ground diffusivity equal to 10−6 m2/s and a BHE diameter of 0.15 m, the simulated time interval exceeds 713 years. The first simulation yields the isoflux g-function, gq. This g-function is used as time-dependent isothermal BC at the BHE-ground interface for the second simulation. Then, the dimensionless heat flux per unit length at the BHE-ground interface, averaged along the BHE field and from the initial time instant to the dimensionless time τ*, is determined. The inverse of this heat flux is used as a correction factor of gq, to determine the first-trial isothermal g-function, gT1. Then, gT1 is employed as isothermal BC at the BHE-ground interface in the third simulation, to determine a new time evolution of the length-averaged dimensionless heat flux per unit length at the BHE-ground interface, which yields the correction factor for gT1, and thus the second-trial isothermal g-function, gT2. The iterations can be stopped when the last applied correction factor is sufficiently close to 1, with respect to the desired accuracy. Finally, one can check the accuracy of the isothermal g-function obtained by performing an additional simulation with this g-function as boundary condition, and comparing the length-averaged time dependent dimensionless heat flux per unit length at the BHE surface to the desired value, 1. The numerical simulations are performed by setting the parameters absolute tolerance and relative tolerance in COMSOL Multiphysics both equal to 10−4, in order to minimize the error due to the time

Fig. 3. Plots of grh, gh and gees versus the logarithm of the dimensionless time, field with buried depth 4 m and borehole spacing 7.5 m. 281

Geothermics 77 (2019) 278–287

C. Naldi, E. Zanchini

Table 1 Comparison of grh, gh and gees.

Table 2 Comparison of the isoflux g-functions gq (adiabatic external surface) and gq

log10(τ*)

grh

gh

gees

100 × (grh – gees)/ gees

100 × (gh – gees)/ gees

0 1 2 3 4 5 6

0.1991 0.3572 0.5321 0.8013 1.5016 2.2314 2.4921

0.2025 0.3622 0.5393 0.8122 1.5217 2.2613 2.5255

0.2030 0.3629 0.5403 0.8137 1.5246 2.2656 2.5303

−1.90 −1.56 −1.53 −1.52 −1.51 −1.51 −1.51

−0.24 −0.20 −0.19 −0.19 −0.19 −0.19 −0.19

(isothermal external surface), field with buried depth 4 m and borehole spacing 7.5 m.

extrapolated estimate of the exact solution in Fig. 4. Mesh 2 is adopted since it yields more precise results with a lower computation time. The relative discrepancy between g2 and gees is a decreasing function of time, ranging from –0.31%, for log10 (τ*) = 0, to –0.17% for log10 (τ*) = 6. In the following, g2 will be considered as the resulting isoflux g-function and will be denoted by gq. In order to check if the results are independent of the domain size, the simulation for gq is repeated by replacing the adiabatic BC at the cylindrical outer surface and at the bottom surface of the domain with the isothermal BC, T* = 0. In fact, the adiabatic BC overestimates the gfunction, whereas the isothermal BC underestimates it. In Table 2, the gfunction obtained with the adiabatic BC, gq, is compared with the gfunction obtained with the zero-temperature BC, gq , for several values of the logarithm of the dimensionless time. The maximum relative discrepancy between the two g-functions is 0.06% and occurs at the last time instant. This result clearly shows that the domain extension is adequate.

log10 (τ*)

gq

gq

100 × (gq

0 1 2 3 4 5 6

0.2023 0.3619 0.5390 0.8119 1.5217 2.2615 2.5259

0.2023 0.3619 0.5390 0.8119 1.5217 2.2615 2.5243

0.00 0.00 0.00 0.00 0.00 0.00 0.06

gq)/gq

isothermal g-function gT5. The figure shows that q¯l* becomes closer and closer to the desired constant value 1. Plots of the isothermal g-functions gT1, gT2, gT3, gT4, gT5 versus the logarithm of the dimensionless time are illustrated in Fig. 6, in the range 0 ≤ log10 (τ*) ≤ 6. The g-functions are graphically indistinguishable for log10 (τ*) < 4. Slight differences become visible for the highest values of τ*, and are highlighted by the zoom on the graph, which refers to the time interval 5.7 ≤ log10 (τ*) ≤ 6. To have a reference in real time scale, one can consider that, for D = 15 cm and αg = 10−6 m2/s, log10 (τ*) = 4 corresponds to 7.135 years and log10 (τ*) = 6 corresponds to 713.5 years. In Table 3, the values of the isothermal g-functions obtained in each iteration are reported for several values of the logarithm of the dimensionless time, together with the percent difference between the values of the last two g-functions. The table shows that this percent difference is extremely small and does not exceed 0.03% for the selected time instants. The maximum percent difference is 0.046% and occurs at log10 (τ*) = 5.65. Moreover, the maximum absolute difference between the values of the last two g-functions, which occurs at log10 (τ*) = 5.70, is 0.0011. Therefore, further iterations are useless, because they would yield corrections much smaller than the discrepancy between the results obtained with the selected mesh and the extrapolated estimate of the exact solution. It is interesting to note that even the difference between the first and the last g-function is rather small: the maximum absolute difference is 0.023 and occurs at log10 (τ*) = 5.15, the maximum relative difference is 1.04% and occurs at log10 (τ*) = 5.05. Therefore, even the first g-function, obtained by two simulations, can be considered as sufficiently accurate for technical purposes. The last isothermal g-function, gT5, that in the following we will call simply gT, is validated by comparison both with the isothermal g-function evaluated by Cimmino and Bernier (2014) for the same BHE field,

3. Results for the BHE field with buried depth 4 m and BHE spacing 7.5 m The desired boundary condition for the isothermal g-functions is: uniform temperature of the BHE-field surface and dimensionless heat flux per unit area, averaged on the BHE-field surface, equal to 1/π, i.e. length-averaged dimensionless linear heat flux q¯l* = 1. A uniform timedependent temperature of the BHE-field surface is imposed as a boundary condition, and the corresponding time evolution of q¯l* is determined. The convergence of q¯l* to the constant value 1 for the field with buried depth 4 m and BHE spacing 7.5 m is illustrated in Fig. 5, where plots of q¯l* versus the dimensionless time are reported for the isothermal g-functions obtained with the recursive procedure illustrated in Section 2. The first plot, denoted by q¯l1* , refers to the third simulation and shows the dimensionless length-averaged linear heat flux produced by the first isothermal g-function, gT1. The last, denoted by q¯l5*, shows the dimensionless length-averaged linear heat flux produced by the

Fig. 5. Plots of q¯l* versus dimensionless time, for the g-functions gT1, gT2, gT3, gT4, gT5, field with buried depth 4 m and borehole spacing 7.5 m.

Fig. 4. Plots of gees, g1 and g2 versus the logarithm of the dimensionless time, field with buried depth 4 m and borehole spacing 7.5 m. 282

Geothermics 77 (2019) 278–287

C. Naldi, E. Zanchini

Fig. 6. Isothermal g-functions gT1, gT2, gT3, gT4, gT5 versus log10 (τ*), field with buried depth 4 m and borehole spacing 7.5 m.

Fig. 7. Isothermal g-function, gT, compared with that determined by Cimmino and Bernier (2014), gCB, field with buried depth 4 m and borehole spacing 7.5 m.

and with the isothermal g-function obtained by applying the method developed by Monzó et al. (2015), that will be called High-ConductiveMaterial (HCM) method. Our g-function, gT, and that determined by Cimmino and Bernier (2014) with each borehole divided into 12 segments, gCB, are plotted versus log10 (τ*) in Fig. 7. The two g-functions are nearly coincident for most of the time, and show small discrepancies for very low and very high values of the dimensionless time. The small differences for low values of τ* are due to the different schemes employed for the simulations: the FLS scheme in (Cimmino and Bernier, 2014), and the FCS scheme in the present paper. The latter yields higher values of the gfunction, for low values of τ*, as is shown, for instance, by Zanchini and Lazzari (2013). The relative difference between gT and gCB at the last time instant is equal to 1.13%. Cimmino and Bernier (2014) evaluated also a more precise g-function of the field, by dividing the boreholes into 256 segments, and reported in their Table 2 the value of the g-function at 20 years. The difference between this value and the corresponding value obtained by means of the numerical method presented here is 0.32%. The HCM model by Monzó et al. (2015) is implemented in our dimensionless computational domain by filling the BHEs with a fictitious HCM, and by linking the BHEs to a horizontal bar of the same HCM through vertical cylindrical HCM connections. The height of the bar is 1 m (i.e. 6.667 diameters) and the depth is 0.15 m, i.e. 1 diameter. The bar starts in correspondence of the symmetry plane that cuts the central BHE, and ends with a solid half cylinder that merges with the vertical cylinder extending from the lateral BHE, as is illustrated in Fig. 8. The adiabatic boundary condition is set at the outer surfaces of the HCM solid placed above the BHEs, except at the upper surface of the horizontal bar, where the time-constant dimensionless total heat flux is applied, and at the upper surface of the BHEs, where the material continuity is imposed. The mesh adopted in the ground is very similar to that employed in our method, and has 1 382 093 tetrahedral elements. A mesh having

141 834 tetrahedral elements is applied in the HCM, so that the mesh has in total 1 523 297 tetrahedral elements. After some simulations with different pairs of values of the thermal conductivity and of the heat capacity per unit volume of the HCM, we have selected the pair with thermal conductivity equal to 107 times that of the ground and heat capacity per unit volume equal to 10−7 times that of the ground. The comparison between the isothermal g-function gT obtained by our model and that obtained by the HCM model, denoted by gHCM, is illustrated in Fig. 9. The figure shows that gT and gHCM are graphically indistinguishable if the whole dimensionless-time interval between log10 (τ*) = 0 and log10 (τ*) = 6 is represented. Very small differences are highlighted by the zoom on the graph, in the dimensionless-time interval 5.7 ≤ log10 (τ*) ≤ 6. The highest relative discrepancy is 0.25%, and occurs at log10 (τ*) = 0.1. A comparison between the isothermal g-function, gT, and the isoflux g-function, gq, is illustrated in Fig. 10. The two g-functions are graphically indistinguishable for log10 (τ*) < 3.5, i.e., for τ < 2.256 years, if D = 15 cm and αg = 10−6 m2/s. At higher times they diverge, with higher values for gq. The maximum relative difference is about 6.6% and occurs at the last time instant, i.e. for log10 (τ*) = 6. The result confirms that the adoption of an isoflux BC at the surface of the BHE field, instead of an isothermal BC, can yield a significant overestimation of the g-function for high values of time. The boundary condition adopted for the isothermal g-functions yields a non-uniform distribution of the heat flux among the boreholes of the field. Fig. 11 illustrates the time evolution of the length-averaged

Table 3 Values of gT1, gT2, gT3, gT4, gT5 at selected time instants, and percent difference between gT5 and gT4, field with buried depth 4 m and borehole spacing 7.5 m. log10 (τ*)

gT1

gT2

gT3

gT4

gT5

100 × (gT5 – gT4)/ gT5

0 1 2 3 4 5 6

0.2020 0.3615 0.5381 0.8099 1.5018 2.1748 2.3794

0.2021 0.3615 0.5381 0.8096 1.4959 2.1555 2.3630

0.2020 0.3615 0.5381 0.8094 1.4943 2.1522 2.3659

0.2021 0.3615 0.5381 0.8094 1.4940 2.1521 2.3685

0.2021 0.3615 0.5381 0.8094 1.4940 2.1524 2.3693

0.00 0.00 0.00 0.00 0.00 0.01 0.03

Fig. 8. Connection between the HCM horizontal bar and the HCM vertical cylinder extending from the lateral BHE, with mesh, field with buried depth 4 m and borehole spacing 7.5 m. 283

Geothermics 77 (2019) 278–287

C. Naldi, E. Zanchini

Fig. 12. Dimensionless temperature averaged on the circumference of the central boreholes, versus z*, for isoflux and isothermal BC, field with buried depth 4 m and borehole spacing 7.5 m.

Fig. 9. Isothermal g-function, gT, compared with that obtained by applying the HCM method, gHCM, field with buried depth 4 m and borehole spacing 7.5 m.

Fig. 12 illustrates the distribution along the BHE length of the dimensionless surface temperature, T*, averaged on the circumference of the central BHEs, for the isoflux BC and for the isothermal BC, at selected time instants. For log10 (τ*) = 0, T* given by the isoflux BC (black dotted line in Fig. 12) is almost uniform, with the same value of T*as that of the isothermal BC (orange dotted line). At higher times (log10 (τ*) = 3), the distribution of T* yielded by the isoflux BC is still uniform, except close to the BHE top (z* = 0) and bottom (z* = 1000), where T* decreases (see the black dashed line in Fig. 12). At the last time instant (log10 (τ*) = 6), the difference between the two temperature distributions is large (see black and orange solid lines in Fig. 12). The isoflux BC yields an inverted-U profile of T*, with values at the BHE top lower than those at the bottom. Fig. 13 shows the distribution along z* of the dimensionless linear heat flux ql* on the surface of the central BHEs, for the isothermal BC, at selected time instants. At log10 (τ*) = 0, the heat flux is nearly uniform, with small peaks at the BHE ends (see dotted line in Fig. 13). The heat flux distribution at log10 (τ*) = 3 is similar to that at log10 (τ*) = 0, with ql* slightly lower for intermediate values of z* and lateral peaks affecting a slightly wider range of z*. At the last time instant (log10 (τ*) = 6), ql* has a definite U-profile, with a peak at the BHE top (z* = 0) higher than that at the bottom (z* = 1000).

Fig. 10. Isoflux g-function, gq, and isothermal g-function, gT, versus log10 (τ*), field with buried depth 4 m and borehole spacing 7.5 m.

4. Results for the BHE field with buried depth 1.5 m and BHE spacing 7.5 m In order to study the influence of the buried depth on the difference

Fig. 11. Length-averaged dimensionless linear heat flux for the central boreholes, q¯la* , and for the lateral boreholes, q¯lb* , versus log10 (τ*), field with buried depth 4 m and borehole spacing 7.5 m.

dimensionless linear heat flux for the central boreholes of the 3 × 2 field, q¯la* , and for the lateral boreholes, q¯lb*. For log10 (τ*) < 2.35, q¯la* and q¯lb* are nearly constant and equal to 1. For log10 (τ*) > 2.35, the lengthaveraged heat flux of the lateral BHEs increases with time, whereas that of the central BHEs decreases, since interference effects become significant and are more important for the central BHEs. For values of log10 (τ*) higher than 5, the length-averaged dimensionless linear heat fluxes are constant in time, and q¯lb* is about 1.25 times q¯la* .

Fig. 13. Plots of ql* versus z* for the central boreholes, with isothermal BC, field with buried depth 4 m and borehole spacing 7.5 m. 284

Geothermics 77 (2019) 278–287

C. Naldi, E. Zanchini

Fig. 16. Isothermal g-functions of the field with buried depth 4 m and of the field with buried depth 1.5 m, BHE spacing 7.5 m.

Fig. 14. Plots of q¯l* versus dimensionless time, for the g-functions gT1, gT2, gT3, gT4, gT5, field with buried depth 1.5 m and BHE spacing 7.5 m.

The maximum absolute difference between the first and the last gfunction is 0.029 and occurs at log10 (τ*) = 5.10; the maximum relative difference is 1.36% and occurs at log10 (τ*) = 5.00. This confirms that even the first g-function can be sufficiently accurate for technical purposes. A comparison between the isothermal g-function obtained by the last iteration, gT, and the isoflux g-function, gq, is illustrated in Fig. 15. As in the case with buried depth 4 m, the two g-functions are graphically indistinguishable for log10 (τ*) < 3.5, whereas diverge at higher times, with higher values for gq. The relative difference is about 6.73% for log10 (τ*) = 5 (i.e. τ = 71.35 years for D = 15 cm and αg = 10−6 m2/s). The maximum relative difference is 8.67% and occurs at the last time instant. The difference between gq and gT is always higher than in the case with Bd = 4 m. A decrease of the buried depth yields a decrease of both the isoflux and the isothermal g-function. The decrease of the isoflux g-function is small and reaches its maximum value, 1.4%, at the last time instant. The effect of the decreased buried depth on the isothermal g-function is more significant and is evidenced in Fig. 16, where the isothermal gfunctions obtained with the different values of Bd are compared. The decrease of the isothermal g-function is graphically visible for log10 (τ*) > 3.5 and reaches its maximum value, 3.26%, at the last time instant.

Table 4 Values of gT1, gT2, gT3, gT4, gT5 at selected time instants, and percent difference between gT5 and gT4, field with buried depth 1.5 m and BHE spacing 7.5 m. log10 (τ*)

gT1

gT2

gT3

gT4

gT5

100 × (gT5 – gT4)/ gT5

0 1 2 3 4 5 6

0.2020 0.3615 0.5381 0.8086 1.4866 2.1241 2.3041

0.2021 0.3615 0.5380 0.8081 1.4788 2.0991 2.2839

0.2020 0.3615 0.5380 0.8079 1.4764 2.0950 2.2878

0.2021 0.3615 0.5380 0.8079 1.4759 2.0951 2.2911

0.2021 0.3615 0.5380 0.8078 1.4759 2.0956 2.2920

0.00 0.00 0.00 −0.01 0.00 0.02 0.04

between isothermal and isoflux g-functions, simulations are performed, for a 3 × 2 BHE field with the same geometrical parameters as those considered in Section 3, except for Bd = 1.5 m (B*d = 10). This is a typical value of the buried depth for Italy, whereas the value Bd = 4 m is typical for countries with colder climates. The recursive procedure described in Section 2 to evaluate isothermal g-functions is applied, and the convergence of the length-averaged dimensionless linear heat flux on the BHE-field surface, q¯l*, is shown in Fig. 14. The figure shows that, also in this case, q¯l* gets closer and closer to the desired constant value 1. The values at selected time instants of the isothermal g-functions obtained in each iteration are reported in Table 4, together with the percent difference between the values of the last two g-functions. The maximum relative difference is 0.057% and occurs at log10 (τ*) = 5.60.

5. Results for the BHE field with buried depth 1.5 m and BHE spacing 6 m To analyze the effect of the borehole spacing on the difference between isothermal and isoflux g-functions, simulations are performed, for a 3 × 2 BHE field with the same geometrical parameters as those considered in Section 4 (D = 0.15 m, H = 150 m, Bd = 1.5 m), except for a borehole spacing S = 6 m. A comparison between the isothermal g-function, gT, and the isoflux g-function, gq, for this borehole field, is illustrated in Fig. 17. The relative difference between gq and gT is nearly the same as that occurring for the BHE field examined in Section 4. A decrease of the BHE spacing yields an increase of both the isoflux and the isothermal g-function. The relative increase is nearly the same for gq and gT, grows with time and equals about 4% at the last time instant. A comparison between the isothermal g-functions of the three BHE fields considered is illustrated in Fig. 18. The figure shows that the values of the isothermal g-function increase when either the buried depth increases or the BHE spacing decreases. The effect of reducing the BHE spacing is more relevant than that of increasing the buried depth.

Fig. 15. Isoflux g-function, gq, and isothermal g-function, gT, versus log10 (τ*), field with buried depth 1.5 m and BHE spacing 7.5 m. 285

Geothermics 77 (2019) 278–287

C. Naldi, E. Zanchini

constant wall heat flux yields an overestimation of the g-function increasing with time. The maximum relative difference found between gq and gT is 6.6% for the field with buried depth 4 m and BHE spacing 7.5 m, and about 8.7% for the other fields. It occurs at the last dimensionless-time instant considered, that corresponds to about 713 years for a borehole diameter of 15 cm and a ground thermal diffusivity of 10−6 m2/s. Both an increase of the buried depth and a decrease of the BHE spacing yield an increase of the isothermal g-function. The effect of reducing the BHE spacing is more relevant than that of increasing the buried depth. Acknowledgements We would like to thank Prof. M. Cimmino for sending us the numerical values of the g-function for the 3 × 2 BHE field plotted in Cimmino and Bernier (2014).

Fig. 17. Isoflux g-function, gq, and isothermal g-function, gT, versus log10 (τ*), field with buried depth 1.5 m and BHE spacing 6 m.

References ASHRAE, 2007. Handbook, HVAC Applications. Chapter 32.. . Bandos, T.V., Montero, A., Fernandez, E., Santander, J.L.G., Isidro, J.M., Perez, J., Fernandez de Cordoba, P.J., Urchueguía, J.F., 2009. Finite line-source model for borehole heat exchangers: effect of vertical temperature variations. Geothermics 38, 263–270. Cimmino, M., 2015. The effects of borehole thermal resistances and fluid flow rate on the g-functions of geothermal bore fields. Int. J. Heat Mass Transfer 91, 1119–1127. Cimmino, M., 2018. Fast calculation of the g-functions of geothermal borehole fields using similarities in the evaluation of the finite line source solution. J. Build. Perform. Simul. https://doi.org/10.1080/19401493.2017.1423390. Cimmino, M., Bernier, M., 2014. A semi-analytical method to generate g-functions for geothermal bore fields. Int. J. Heat Mass Transfer 70, 641–650. Cimmino, M., Bernier, M., Adams, F., 2013. A contribution towards the determination of g-functions using the finite line source. Appl. Therm. Eng. 51, 401–412. Claesson, J., Eskilson, P., 1987. Conductive Heat Extraction by a Deep Borehole, Analytical Studies. Technical Report. Lund University, Sweden. Claesson, J., Javed, S., 2011. An analytical method to calculate borehole fluid temperatures for time-scales from minutes to decades. ASHRAE Trans. 117 (2), 279–288. Costes, V., Peysson, P., 2008. Capteurs géothermiques verticaux enterrés: Validation expérimentale de nouveaux modèles développés dans l’environnement TRNSYS. École Polytechnique de Montréal, Montréal, pp. 89. Eskilson, P., 1987. Thermal Analysis of Heat Extraction Boreholes. Ph.D. Thesis. University of Lund, Lund, Sweden. Fisher, D.E., Rees, S.J., Padhmanabhanc, S.K., Murugappan, A., 2006. Implementation and validation of ground-source heat pump system models in an integrated building and system simulation environment. HVAC&R Res. 12, 693–710. Fossa, M., 2011. The temperature penalty approach to the design of borehole heat exchangers for heat pump applications. Energy Build. 43, 1473–1479. Fossa, M., Cauret, O., Bernier, M., 2009. Comparing the thermal performance of ground heat exchangers of various lengths. Proceedings of Effstock 2009. Stockholm, Sweden, pp. 8. Hellström, G., Sanner, B., 1994. Software for dimensioning of deep boreholes for heat extraction. Proceedings of Calorstock 1994. Espoo/Helsinki, Finland, pp. 195–202. Kavanaugh, S.P., Rafferty, K., 1997. Ground-Source Heat Pumps. Design of Geothermal Systems for Commercial and Institutional Buildings. ASHRAE, Atlanta, GA, USA. Lamarche, L., 2017. G-function generation using a piecewise-linear profile applied to ground heat exchangers. Int. J. Heat Mass Transfer 115, 354–360. Lamarche, L., Beauchamp, B., 2007. A new contribution to the finite line-source model for geothermal boreholes. Energy Build. 39, 188–198. Liu, X., Hellström, G., 2006. Enhancements of an integrated simulation tool for groundsource heat pump system design and energy analysis. Proc. 10th International Conference on Thermal Energy Storage. Richard Stockton College of New Jersey. Lund, J.W., Boyd, T.L., 2016. Direct utilization of geothermal energy 2015 worldwide review. Geothermics 60, 66–93. Monzó, P., Acuña, J., Fossa, M., Palm, B., 2013. Numerical Generation of the Temperature Response Factors for a Borehole Heat Exchangers Field. European Geothermal Congress 2013, Pisa, Italy, pp. 1–8. Monzó, P., Mogensen, P., Acuña, J., Ruiz-Calvo, F., Montagud, C., 2015. A novel numerical approach for imposing a temperature boundary condition at the borehole wall in borehole fields. Geothermics 56, 35–44. Monzó, P., Puttige, A.R., Acuña, J., Mogensen, P., Cazorla, A., Rodriguez, J., Montagud, C., Cerdeira, F., 2018. Numerical modeling of ground thermal response with borehole heat exchangers connected in parallel. Energy Build. 172, 371–384. Naldi, C., Zanchini, E., 2017. Hourly simulation of a Ground-coupled heat pump system. J. Phys. Conf. Ser. 796, 10 012020. Naldi, C., Zanchini, E., 2018. Effects of the total borehole length and of the heat pump inverter on the performance of a ground-coupled heat pump system. Appl. Therm. Eng. 128, 306–319. Oberkampf, W.L., Roy, C.J., 2010. Verification and Validation in Scientific Computing. Cambridge University Press, New York.

Fig. 18. Isothermal g-functions of the 3 × 2 BHE fields considered.

6. Conclusions We have presented a new numerical method to determine isothermal g-functions of borehole heat exchanger (BHE) fields, i.e. gfunctions with the boundary condition of uniform temperature and timeconstant heat flux per unit length averaged along BHE field. The method is iterative. The g-function gq obtained by the boundary condition of uniform and constant heat flux is used as a boundary condition to determine the first-trial time evolution of the dimensionless lengthaveraged linear heat flux q¯l1*. The inverse of the time average of q¯l1* from the initial instant to the current time instant is used as a correction factor of gq, to obtain the first-trial g-function with uniform wall temperature, gT1. Then, gT1 is used as a boundary condition to determine the length-averaged linear heat flux that yields the second-trial gfunction with uniform wall temperature, gT2, and so on. The first-trial gfunction is already accurate enough for technical purposes, while several iterations can yield very accurate benchmark isothermal g-functions. The method has been applied to compare the isoflux and the isothermal g-function for three 3 × 2 borehole fields, with the same BHE length (150 m) and diameter (0.15 m), and with the following values of the other geometrical parameters: buried depth 4 m and BHE spacing 7.5 m, buried depth 1.5 m and BHE spacing 7.5 m, buried depth 1.5 m and BHE spacing 6 m. The isothermal g-function for the field with buried depth 4 m and BHE spacing 7.5 m has been validated by comparison with that determined by Cimmino and Bernier (2014) and with that obtained by applying the method developed by Monzó et al. (2015). The results confirm that the boundary condition of uniform and 286

Geothermics 77 (2019) 278–287

C. Naldi, E. Zanchini Priarone, A., Fossa, M., 2016. Temperature response factors at different boundary conditions for modelling the single borehole heat exchanger. Appl. Therm. Eng. 103, 934–944. Puttige, A.R., Rodriguez, J., Monzó, P., Cerdeira, F., Fernández, A., Novelle, L., Acuña, J., Mogensen, P., 2016. Improvements on a Numerical Model of Borehole Heat Exchangers. European Geothermal Congress 2016, Strasbourg, France, pp. 10. Richardson, L.F., 1911. The approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stresses in a masonry dam. Philos. Trans. R. Soc. Lond. Ser. A Containing Papers Math. Phys. Charact. 210, 307–357. Richardson, L.F., 1927. The deferred approach to the limit. Part I. Single lattice. Philos. Trans. R. Soc. Lond. Ser. A Containing Papers Math. Phys. Charact. 226, 299–349.

Self, S.J., Reddy, B.V., Rosen, M.A., 2013. Geothermal heat pump systems: Status review and comparison with other heating options. Appl. Energy 101, 341–348. Spitler, J.D., 2000. GLHEPRO – a design tool for commercial building ground loop heat exchanger. Proceedings of the 4th International Conference on Heat Pumps in Cold Climates. Aylmer, QC, Canada, pp. 1–16. Zanchini, E., Lazzari, S., 2013. Temperature distribution in a field of long borehole heat exchangers (BHEs) subjected to a monthly averaged heat flux. Energy 59, 570–580. Zanchini, E., Lazzari, S., 2014. New g-functions for the hourly simulation of double U-tube borehole heat exchanger fields. Energy 70, 444–455. Zeng, H.Y., Diao, N.R., Fang, Z.H., 2002. A finite line-source model for boreholes in geothermal heat exchangers. Heat Transfer–Asian Res. 31, 558–567.

287