New g-functions for the hourly simulation of double U-tube borehole heat exchanger fields

Energy xxx (2014) 1e12

Contents lists available at ScienceDirect

Energy journal homepage: www.elsevier.com/locate/energy

New g-functions for the hourly simulation of double U-tube borehole heat exchanger fields E. Zanchini*, S. Lazzari University of Bologna, Department of Industrial Engineering (DIN), Bologna, Italy

a r t i c l e i n f o

a b s t r a c t

Article history: Received 23 January 2014 Received in revised form 13 March 2014 Accepted 6 April 2014 Available online xxx

A new method for the hourly simulation of fields of long BHEs (borehole heat exchangers) is presented, in dimensionless form. The method considers also the internal structure of the BHE, with reference to double U-tube BHEs. The analytical expressions of new g-functions (time evolutions of the dimensionless temperature, averaged along the BHE length, due to a uniform and constant heat load) are presented, which complement those reported in Zanchini E., Lazzari S., Energy 2013; 59: 570e80. The new gfunctions yield the dimensionless temperature at the interface tubes-grout, and allow a more precise determination of the time evolution of the temperature of the operating fluid during hourly peak loads. Moreover, g-functions at the interface tubes-grout and at several radial distances are reported also for BHEs with a ratio 500 between length and diameter, a value not considered in the previous paper. Due to its higher accuracy, the new method is recommended for the hourly simulation and the optimization of ground-coupled heap pump systems with double U-tube BHEs. Ó 2014 Elsevier Ltd. All rights reserved.

Keywords: Ground-coupled heat pumps BHE (borehole heat exchanger) fields Hourly simulation g-functions Finite element analysis

1. Introduction GSHPs (ground-source heat pumps) are a very efficient technology for building heating and cooling. If the heat exchanger with the ground is a closed-loop one, GSHPs are called GCHPs (groundcoupled heat pumps) [1]. The most widely used heat exchangers for GCHPs are vertical, and are called BHEs (borehole heat exchangers); they are usually composed of either a single U-tube or a double Utube, or two coaxial tubes. In this paper, double U-tube BHEs are considered, which are the most widely employed. A wide research work on the thermal performance and sustainability of BHE fields is available in the literature. Relevant contributions concern the evaluation of TRTs (thermal response tests) [2e7], the analysis of the internal heat transfer mechanisms in a single BHE [8e12], the long-term sustainability of BHE fields in the absence of groundwater movement [13,14], the effects of the groundwater movement [15e17]. Some recent papers analyze the thermal response of pile heat exchangers, where heat transfer pipes are cast into the building piled foundations [18e20]. Here, we focus our attention on the available design methods of BHE fields. Pioneering contributions to the design of BHE fields were given by Claesson and Eskilson [21,22], who introduced the concept of g-

* Corresponding author. E-mail address: [email protected] (E. Zanchini).

function, i.e. the dimensionless thermal response due to a constant heat load, and suggested the use of a semi-analytical expression of the temperature field produced by a FLS (finite line source) subjected to a constant heat flux per unit length. A simple and widely employed design method for BHE fields was developed by Kavanaugh and Rafferty [23] and is recommended by ASHRAE [1]. It is based on the solution of the equation for the heat transfer from an infinitely long cylinder placed in a homogeneous solid medium, and considers the superposition of three heat pulses, each with a constant power, which account respectively for seasonal heat imbalances, average monthly heat load during the design month, and peak heat pulse during the design day. The method takes into account the thermal interference between BHEs, but is limited to 10 years of operation and, therefore, does not guarantee the longterm sustainability for several decades. A commercial design code for BHE fields which grants the longterm sustainability, called EED (earth energy designer) was developed by Hellström and Sanner [24]. EED is based on expressions of the dimensionless temperature produced by several configurations of BHE fields (g-functions), obtained through numerical simulations. Zeng, Diao and Fang [25] pointed out that the use of the semianalytical expression of the temperature field produced by a FLS suggested by Claesson and Eskilson, evaluated at the middle of the length of the BHE, yields an overestimation of the mean temperature field at the BHE surface; they proposed to employ the value given by that expression when averaged along the BHE length,

http://dx.doi.org/10.1016/j.energy.2014.04.022 0360-5442/Ó 2014 Elsevier Ltd. All rights reserved.

Please cite this article in press as: Zanchini E, Lazzari S, New g-functions for the hourly simulation of double U-tube borehole heat exchanger fields, Energy (2014), http://dx.doi.org/10.1016/j.energy.2014.04.022

2

E. Zanchini, S. Lazzari / Energy xxx (2014) 1e12

Nomenclature

Rt

Ai

S S* T T* Tf Tg Ts u x x 0, x 1 Y0, Y1

dimensionless heat load per unit length during the i-th hour a0, ., a6 dimensionless coefficients b0, ., b6 dimensionless coefficients d distance between opposite tubes (m) D diameter of the BHE (m) Dt external diameter of each tube (m) d* d/D, dimensionless distance between opposite tubes D*t Dt/D, dimensionless external diameter of each tube g g-function H heaviside unit-step function J0, J1 Bessel functions of the first kind, of order 0 and 1 kg thermal conductivity of the ground (W m1 K1) kgt thermal conductivity of the grout (W m1 K1) k* kgt/kg, dimensionless thermal conductivity of the grout L BHE length (m) L* L/D, dimensionless BHE length n outward unit normal Q heat flux per unit length (W m1) Q0 reference heat flux per unit length (W m1) Q* Q/Q0, dimensionless heat flux per unit length r radial coordinate (m) r* r/D, dimensionless radial coordinate rd radius of the computational domain (m) Rb thermal resistance per unit length of the BHE (m K W1) Rgt thermal resistance per unit length of the grout (m K W1)

which has the form of a double integral. Lamarche and Beauchamp [26], by analytical methods, determined a simpler form of that solution, which requires a lower computational time. Bandos et al. [27] determined another expression of the same solution, in the form of a single integral, which, through numerical integrations, yields very accurate results. Moreover, they analyzed the effects of the surface temperature oscillations and of the geothermal gradient. By employing the results obtained in Ref. [27], Fossa [28,29] proposed simple approximate expressions for the temperature averaged along the BHE length, with empirical coefficients determined through the analysis of several different BHE fields. Bernier et al. [30] studied, in dimensionless form, the temperature penalty Tp, which represents an effective increase or decrease of the undisturbed ground temperature to account for BHE thermal interactions. They reported dimensionless tables of Tp, for several BHE field geometries, and an approximate correlation, valid for any BHE field up to a 12
12 square field. Cimmino et al. [31] pointed out that, while in the customary design methods based on the FLS solution the BHEs of a field are considered as subjected to the same power per unit length, in practice they are fed in parallel and have the same mean temperature. As a consequence, the heat extraction (or injection) rate of the central BHEs can differ from that of the peripheral ones. The authors proposed a new method, based on the FLS model, which takes into account the variation of the heat extraction rate among boreholes of the same bore field while keeping the same mean borehole wall temperature for every borehole in the field. Bernier et al. [32] presented a MLAA (multiple load aggregation algorithm) to aggregate heating/cooling loads when using the cylindrical heat source model to perform hourly energy simulations

thermal resistance per unit length of the tubes (m K W 1) surface (m2) surface in dimensionless coordinates temperature ( C) kg(T  Tg)/Q0, dimensionless temperature mean temperature of the fluid ( C) undisturbed ground temperature ( C) mean temperature of the BHE external surface ( C) dimensionless integration variable log10s* threshold values of x Bessel functions of the second kind, of order 0 and 1

Greek symbols V Nabla operator (m1) V* D V, dimensionless Nabla operator (r c)g heat capacity per unit volume of the ground (J m3) (r c)gt heat capacity per unit volume of the grout (J m3) (r c)* (r c)gt/(r c)g, dimensionless heat capacity of the grout s time (s) s* kgs/((r c)g D2), dimensionless time * s0 dimensionless duration of 1 h Subscripts 0.5 at the external surface of the BHE 2D obtained by a 2-D simulation on a cross section FCS obtained by the finite-cylindrical-source model t-gt at the interface tubes-grout

of GCHP systems. They compared the results of their method with those obtained by means of the DST (duct storage) model, implemented in TRNSYS, and found a good agreement. Fossa and Minchio [33] performed a similar comparison, by implementing a modified MLAA method into a model able to employ suitable g-functions generated starting from the FLS solution. They found a good agreement between their method and the DST model for compact BHE fields, which better fulfill the assumptions of the DST model, and some discrepancies for single-line fields. In a recent paper [34], accurate analytical expressions of the gfunctions, in the form of polynomial functions of the logarithm of the dimensionless time, were determined for fields of BHEs with length L and diameter D, for the following values of the ratio L/D: 2000, 1400, 1000, 700. The scheme of a FCS (finite cylindrical source) was used for the BHEs, and, for each value of L/D, coefficients of the polynomial g-functions were reported for several values of the dimensionless radial distance r/D, starting from the BHE surface (r/D ¼ 0.5). All the design methods quoted above, including the tables of gfunctions reported in Ref. [34], cannot be considered as completely satisfactory for accurate hourly simulations of BHE fields. In fact, they allow to determine with a good accuracy the time evolution of the temperature Ts of the external surface of the BHEs. To determine the time evolution of the mean temperature Tf of the working fluid, one should then employ the quasi-stationary approximation for the heat transfer inside each BHE, which relates Tf to Ts by the equation

Tf ¼ Ts þ Rb Q ;

(1)

where Rb is the BHE thermal resistance per unit length and Q is the heat load per unit length, positive if heat is supplied to the ground.

Please cite this article in press as: Zanchini E, Lazzari S, New g-functions for the hourly simulation of double U-tube borehole heat exchanger fields, Energy (2014), http://dx.doi.org/10.1016/j.energy.2014.04.022

E. Zanchini, S. Lazzari / Energy xxx (2014) 1e12

3

However, for a fast varying heat load, the heat transfer in the grout is not quasi-stationary and Eq. (1) can yield non-negligible errors. In order to allow a more precise hourly simulation of BHE fields, even in the presence of fast varying heat loads, in this paper new gfunctions are presented, which yield directly the temperature (averaged in the BHE length) at the interface between the tubes which contain the working fluid and the grout of the BHE, which will be denoted by Ttegt. The evaluations were performed for double U-tube BHEs, with several different values of the properties of grout and ground. If Ttegt is known, the mean temperature of the working fluid can be determined as

Tf ¼ Ttgt þ Rt Q ;

(2)

where Rt is the thermal resistance of the (polyethylene) tubes, including the convective thermal resistance between fluid and tubes. The g-functions for the evaluation of Ttegt, denoted by gtegt, were determined by the following procedure. For low values of the dimensionless time, where the effects of the axial heat transfer are negligible, the new g-functions were determined through 2-D finite element simulations in which the real cross section of the BHE is considered. For high values of the dimensionless time, the BHE was sketched as a finite-length cylindrical heat source and 2-D axisymmetric simulations were employed; the g-function at the interface tubes-grout, gtegt, was obtained by adding to the g-function at the outer surface of the BHE, g0.5, the asymptotic value of the difference gtegt  g0.5, determined through the 2-D simulations of the cross section. Contrary to those at the interface tubes-grout, the g-functions at a distance from the BHE, which are necessary to take into account the interference between BHEs, are practically independent of the internal structure of the BHE. In particular, as is shown in Ref. [34], the g-functions obtained by the FCS (finite-cylindrical-source) model are graphically indistinguishable from those obtained by the FLS model. Therefore, the g-functions reported in Ref. [34], for r/ D ¼ 30, 40, 60, 80, 120, 170, 230, 400, 600 can be used for any internal structure of the BHEs and any time evolution of the heat load. These g-functions were determined for L/D ¼ 2000, 1400, 1000, 700; those for L/D ¼ 500, obtained by the same method (finite-element simulations performed with the FCS model) are added here. The paper is organized as follows. In Section 2, we describe the method employed to determine the g-functions at the interface tubes-grout for low values of the dimensionless time, and the validations of the method. In Section 3, we present the results on the g-functions at the interface tubes-grout for low values of the dimensionless time. In Section 4 we describe the method employed to determine the g-functions at the interface tubes-grout for both low and high values of the dimensionless time. In Section 5 we present the final results on the new g-functions at the interface tubes-grout. In Section 6 we present the results on the g-functions at the external surface of the BHE for L/D ¼ 500, obtained by applying the finite-cylindrical-source model, which complement those reported in Ref. [34]. Finally, in Section 7 we illustrate, by an example, the errors that can be caused by simulation methods based on the g-functions at the external surface of the BHE and the quasi-stationary approximation of the heat transfer in the grout.

Fig. 1. Cross section of the BHE grout, with dimensionless lengths.

dimensionless parameters which determine the cross section of the BHE grout (the most relevant part of the BHE for the effects of thermal inertia), are D*t ¼ Dt =D and d* ¼ d/D. The following values of these parameters, taken from a widely employed double U-tube BHE with diameter D ¼ 152 mm, have been adopted in computations: D*t ¼ 0:2106 (i.e., Dt ¼ 32 mm), d* ¼ 0.5592. A sketch of the BHE cross section considered, with dimensionless lengths, is reported in Fig. 1. Let us define the dimensionless temperature T* and the dimensionless time s* as follows:

T * ¼ kg

T  Tg ; Q0

s* ¼

kg s ; ðrcÞg D2

(3)

where T is temperature, s is time, kg and (rc)g are the thermal conductivity and the heat capacity per unit volume of the ground, Tg is the undisturbed ground temperature, and Q0 is a reference heat load per unit length. In analogy with Ref. [34], we will determine the dimensionless temperature, averaged along the BHE length, due to a constant heat load per unit length with value Q0, which begins at the initial instant s ¼ 0; as usual, this quantity will be called g-function. As is well known, the effects of any stepwise variable heat load supplied to a BHE, at a given position, can be obtained by a superposition of g-functions (at that position) properly displaced in time, and the effects of the other BHEs of the field can be determined by a superposition of g-functions evaluated at different radial distances. As is shown in Figs. 5 and 6 of Ref. [34], for values of the dimensionless time lower than 103 the g-functions can be considered as

2. Evaluation of the g-functions at the interface tubes-grout, for low values of the dimensionless time: method and validations Let us consider a double U-tube BHE with diameter D and length L. Let us denote by Dt the external diameter of each tube and by d the distance between the axes of opposite tubes. The

Fig. 2. Plots of gtegt,2D for k* ¼ 1 and (rc)* ¼ 0.4, 0.7, 1, in the range 4  x  3.

Please cite this article in press as: Zanchini E, Lazzari S, New g-functions for the hourly simulation of double U-tube borehole heat exchanger fields, Energy (2014), http://dx.doi.org/10.1016/j.energy.2014.04.022

4

E. Zanchini, S. Lazzari / Energy xxx (2014) 1e12

Fig. 3. Plots of gtegt,2D  g0.5,2D for k* ¼ 1 and (rc)* ¼ 0.4, 0.7, 1, in the range 4  x  3.

independent of the BHE length L, so that they can be determined by 2-D numerical simulations performed on a cross section of the BHE. Therefore, for low values of s*, we will consider the real internal structure of the BHE and will determine the g-functions at the interface tubes-grout and at the interface BHE-ground through 2-D finite element simulations performed on a cross section, implemented by means of the software COMSOL Multiphysics. Consider a constant heat flux per unit BHE length which begins at the instant of time s ¼ 0 and is expressed by the equation

Q ðsÞ ¼ Q0 HðsÞ

(4)

where H(s) is Heaviside’s unit-step function. The differential equations to be solved, in grout and ground, are, respectively

Fig. 5. Plots of gtegt,2D  g0.5,2D for k* ¼ 0.7 and (rc)* ¼ 0.4, 0.7, 1, in the range 4  x  3.

where Stegt is the surface between tubes and grout and n denotes the outward unit normal. The external boundary of the computational domain, which is a circle with radius rd, is considered as adiabatic. Continuity conditions hold at the interface grout-ground, and the initial condition is T ¼ Tg on the whole computational domain. By introducing, in addition to T* and s*, the dimensionless operator V*¼D V and the dimensionless quantities

k* ¼

kgt ; kg

ðrcÞ* ¼

(5)

vT * k* *2 * ¼ V T ; * vs ðrcÞ*

kg 2 vT ¼ V T; vs ðrcÞg

(6)

vT * ¼ V*2 T * ; vs*


kgt ðVT$nÞ
S

tgt

¼

Q0 HðsÞ; 4pDt

(7)

Fig. 4. Plots of gtegt,2D for k* ¼ 0.7 and (rc)* ¼ 0.4, 0.7, 1, in the range 4  x  3.

ðrcÞg

;

D*t ¼

Dt ; D

(8)

one can rewrite Eqs. (5)e(7) in the following dimensionless form:

kgt 2 vT ¼ V T; vs ðrcÞgt

where kgt and (rc)gt are the thermal conductivity and the heat capacity per unit volume of the grout. The boundary condition at the surface between tubes and grout is

ðrcÞgt



 V* T * $n


S*tgt

(9)

(10)

¼

  1 s* ; H 4pk* D*t

(11)

where S*tgt denotes the surface between tubes and grout in the dimensionless domain. The dimensionless initial condition is T* ¼ 0 in the whole dimensionless computational domain.

Fig. 6. Plots of gtegt,2D for k* ¼ 0.4 and (rc)* ¼ 0.4, 0.7, 1, in the range 4  x  3.

Please cite this article in press as: Zanchini E, Lazzari S, New g-functions for the hourly simulation of double U-tube borehole heat exchanger fields, Energy (2014), http://dx.doi.org/10.1016/j.energy.2014.04.022

E. Zanchini, S. Lazzari / Energy xxx (2014) 1e12 Table 1 Comparison of the results obtained by three different meshes, for k* ¼ 0.7 and (rc)* ¼ 0.7. x

3 2 1 0 1 2 3 4 5

* Ttub

Ts*

Mesh 1

Mesh 2

Mesh 3

Mesh 1

Mesh 2

Mesh 3

0.016281 0.046496 0.153753 0.317382 0.490448 0.671298 0.853832 1.036627 1.219372

0.016135 0.046686 0.153847 0.317337 0.490539 0.671434 0.854049 1.036844 1.219643

0.015876 0.046564 0.154324 0.317755 0.490613 0.671627 0.854371 1.037322 1.220193

8.601E-06 0.004446 0.052847 0.189289 0.360115 0.540914 0.723633 0.906636 1.089592

1.054E-05 0.004405 0.052859 0.189213 0.360140 0.540967 0.723750 0.906737 1.089732

1.624E-05 0.004442 0.052931 0.189547 0.360134 0.541034 0.723900 0.906999 1.090019

The dimensionless Eqs. (9)e(11), with the adiabatic condition at the external boundary, the continuity conditions at the interface grout-ground and initial condition T* ¼ 0, were solved by means of COMSOL Multiphysics, a commercial code based on the finite element method. A computational domain with external dimensionless radius rd/D ¼ 1000 was employed. The dimensionless time interval 104  s*  105 was considered, in the logarithmic scale 4  log10s*  5, and was divided into intervals of log10s* equal to 0.002, with a total of 4500 intervals. The average dimensionless temperatures at the interface tubes-grout and that at the interface grout-ground, called g-functions at those surfaces, were evaluated; they were denoted, respectively, by gtegt,2D and g0.5,2D, to remind that they were evaluated by 2-D simulations on a plane orthogonal to the BHE axis. Calculations were performed with the following values of k* and of (rc)*: k* ¼ 0.4, 0.7, 1.0; (rc)* ¼ 0.4, 0.7, 1.0, in order to cover the typical range of values of the physical properties of grout and ground. To check the mesh independence of results, the g-functions gte * * gt,2D and g0.5,2D, for k ¼ 0.7 and (rc) ¼ 0.7, were evaluated by three different unstructured triangular meshes: Mesh 1, with 7774 triangular elements; Mesh 2, with 21,696 triangular elements; Mesh 3, with 39,904 triangular elements. The mean square deviation between the results obtained through Mesh 1 and those obtained through Mesh 3 is 6.19
105 for gtegt,2D and 2.97
105 for g0.5,2D; that between the results obtained through Mesh 2 and those obtained through Mesh 3 is 4.34
105 for gtegt,2D and 2.21
105 for g0.5,2D. The comparison between the results obtained by Mesh 1, Mesh 2 and Mesh 3, for a few selected values of x ¼ log10s*, is illustrated in Table 1. Mesh 3 was adopted for final computations. In order to check the adequacy of the size of the computational domain adopted, we compared the results obtained with the condition of adiabatic external boundary (which yields an overestimation of gtegt,2D and g0.5,2D) with those obtained with the condition of isothermal external boundary, with T ¼ Tg (which yields an underestimation of gtegt,2D and g0.5,2D). The comparison was Table 2 Comparison of the results obtained by either adiabatic or isothermal external boundary, for k* ¼ 0.7 and (rc)* ¼ 0.7. x

3 2 1 0 1 2 3 4 5

* Ttub

Ts*

Adiabatic

Isothermal

Adiabatic

Isothermal

0.015876 0.046564 0.154324 0.317755 0.490613 0.671627 0.854371 1.037322 1.220193

0.015878 0.046564 0.154322 0.317752 0.490614 0.671627 0.854373 1.037322 1.220185

1.624E-05 0.004442 0.052931 0.189547 0.360134 0.541034 0.723900 0.906999 1.090019

1.619E-05 0.004442 0.052931 0.189544 0.360136 0.541033 0.723902 0.906998 1.090011

5

performed for k* ¼ 0.7 and (rc)* ¼ 0.7, with Mesh 3. The mean square deviation between the results obtained through the isothermal boundary condition and those obtained through the adiabatic boundary condition is 3.06
107 for gtegt,2D and 2.56
107 for g0.5,2D. The check ensures that the extension of the computational domain is more than sufficient. The comparison between the results obtained with the different boundary conditions, for a few selected values of x ¼ log10s*, is illustrated in Table 2. The adiabatic boundary condition was adopted in final computations. The accuracy of our results was checked also by comparing the values of g0.5,2D obtained through our numerical code, with Mesh 3, for very high values of k* and very low values of (rc)*, with those of the dimensionless temperature at the surface of an infinite cylindrical heat source with dimensionless radius 0.5, determined through the analytical solution presented by Carlslaw and Jaeger [35]. With the dimensionless variables T* and s* defined here, and the dimensionless radial coordinate r* ¼ r/D, this solution has the form

  ZN  *  * 2   Y0 2r u J1 ðuÞ  J0 2r * u Y1 ðuÞ 1  e4s u 1 T * r * ; s* ¼ 2 du: p u2 J12 ðuÞ þ Y12 ðuÞ 0

(12) k*

(rc)*

The following pairs of values of and were considered in finite element simulations: k* ¼ 10 and (rc)* ¼ 0.1, k* ¼ 100 and (rc)* ¼ 0.01, k* ¼ 1000 and (rc)* ¼ 0.001. Clearly, when k* / N and (rc)* / 0 the numerical solution should tend to the analytical one. The results of the check are very satisfactory. In the dimensionlesstime interval 104  s*  105, the root mean square deviations between the numerical and the analytical solution are 4.83
104, 1.54
104 and 1.18
104, respectively for k* ¼ 10 and (rc)* ¼ 0.1, k* ¼ 100 and (rc)* ¼ 0.01, k* ¼ 1000 and (rc)* ¼ 0.001. The comparison between the results obtained through the analytical solution and those obtained by the finite element code, for some selected values of x ¼ log10s*, is illustrated in Table 3. 3. Evaluation of the g-functions at the interface tubes-grout, for low values of the dimensionless time: results Plots of gtegt,2D versus x, for k* ¼ 1 and (rc)* ¼ 0.4, 0.7, 1, are reported in Fig. 2. Although computations were performed for 4  x  5, the plots are represented only in the range 4  x  3, because they are indistinguishable from each other for higher values of x. For x > 3, gtegt,2D increases linearly, with the same slope which appears in the figure for 2  x  3. The figure shows that, for k* ¼ 1, the effect of the specific heat capacity of the grout on the dimensionless temperature at the interface tubes-grout is not very important and becomes negligible for x > 1.

Table 3 Comparison between the values of g0.5 obtained through the analytical solution for an infinite cylindrical heat source and those obtained by the finite element simulation code. x

Analytical

k* ¼ 10 (rc)* ¼ 0.1

k* ¼ 100 (rc)* ¼ 0.01

k* ¼ 1000 (rc)* ¼ 0.001

4 3 2 1 0 1 2 3 4 5

0.003560 0.011051 0.033052 0.089861 0.202907 0.363229 0.541942 0.724504 0.907650 1.090870

0.000000 0.004806 0.027476 0.086036 0.201242 0.362609 0.541496 0.724000 0.906943 1.090008

0.000000 0.009672 0.032330 0.089345 0.202504 0.362932 0.541588 0.723958 0.906942 1.089982

0.000000 0.010247 0.032833 0.089689 0.202616 0.363010 0.541543 0.723947 0.907001 1.090020

Please cite this article in press as: Zanchini E, Lazzari S, New g-functions for the hourly simulation of double U-tube borehole heat exchanger fields, Energy (2014), http://dx.doi.org/10.1016/j.energy.2014.04.022

6

E. Zanchini, S. Lazzari / Energy xxx (2014) 1e12

Fig. 7. Plots of gtegt,2D  g0.5,2D for k* ¼ 0.4 and (rc)* ¼ 0.4, 0.7, 1, in the range 4  x  3.

Plots of gtegt,2D  g0.5,2D versus x, for k* ¼ 1 and (rc)* ¼ 0.4, 0.7, 1, are reported in Fig. 3, in the range 4  x  3. The figure shows that the difference between the dimensionless temperature at the interface tubes-grout and that at the interface grout-ground increases and depends on (rc)* only for x < 1. Then, it reaches a constant value, which is independent of (rc)*. Therefore, for x > 1 or close to 1, the quasi-stationary approximation for the heat transfer in the grout can be employed without relevant errors. For typical values of the BHE diameter and of the thermal diffusivity of the ground, such as D ¼ 0.15 m and kg/(rc)g ¼ 106 m2/s, x ¼ 1 corresponds to a time interval of 62.5 h. Plots of gtegt,2D and of gtegt,2D  g0.5,2D versus x, in the range 4  x  3, for k* ¼ 0.7 and (rc)* ¼ 0.4, 0.7, 1 are reported in Figs. 4 and 5 respectively. The main difference, with respect to the cases with k* ¼ 1, is the higher value reached by gtegt,2D  g0.5,2D and, as a consequence, also by gtegt,2D. Plots of gtegt,2D and of gte * gt,2D  g0.5,2D versus x, in the range 4  x  3, for k ¼ 0.4 and (rc)* ¼ 0.4, 0.7, 1 are reported in Figs. 6 and 7 respectively. Clearly, the values reached by gtegt,2D  g0.5,2D and by gtegt,2D are even higher. 4. Evaluation of the g-functions at the interface tubes-grout, for both low and high values of the dimensionless time, 10L4 £ s* £ 106: method As is well known [26,27,34], for high values of the dimensionless time the effects of the heat transfer between the BHE and the ground surface cannot be neglected. Therefore, in order to consider the real internal structure of the BHE, rather heavy 3-D computations should be performed. However, we have shown in Section 3 that, for x > 1 (i.e., s* > 10), the difference between gtegt,2D and g0.5,2D is a constant, which can be determined by assuming steady-state heat transfer in the grout, and is a result of the 2-D simulations already performed. Therefore, for x > 1 (or close to 1) it is possible to determine the time evolution of gtegt by adding to the time evolution of the g-function at the interface BHE-ground, obtained by the FCS model and denoted by g0.5,FCS, the constant difference gtegt,2D  g0.5,2D. In order to determine the best value of x above which the solution for large values of x should replace the 2-D solution, for each set of values of the dimensionless BHE length L* ¼ L/D, of k* and of (rc)*, we employed the following method. We compared gtegt,2D with the g-

Fig. 8. Comparison of gtegt,2D and gtegt,FCS, for L * ¼ 700, k* ¼ 0.7 and (rc)* ¼ 0.7, in the range 4  x  5.

function obtained as a sum of g0.5,FCS and the asymptotic value of gte  g0.5,2D, which will be denoted by gtegt,FCS. For low values of x, gtegt,FCS is greater than gtegt,2D, because the difference between the dimensionless temperature at the interface tubes-grout and that at the interface BHE-ground is overestimated. On the contrary, for high values of x, gtegt,FCS is lower than gtegt,2D, because it takes into account the heat transfer between the BHE and the ground surface. The optimal value for the transition from gtegt,2D to gtegt,FCS is the value of x at which these g-functions become equal. We will denote by gtegt the g-function at the interface tubes-grout obtained by this method, which depends on L*, k*, (rc)*. For L* ¼ 2000, 1400, 1000, 700, the g-functions g0.5,FCS at the interface BHE-ground were already available, since they were determined in Ref. [34]. The results obtained there were employed. In addition, the g-function g0.5,FCS for L* ¼ 500 was determined here, by the same method, namely by finite-element simulations performed through COMSOL Multiphysics. Some details of the simulation code are given in Section 6. An example of application of the method to obtain gtegt is illustrated in Fig. 8, where gtegt,2D and gtegt,FCS are compared, in the range 4  x  5, for L* ¼ 700, k* ¼ 0.7 and (rc)* ¼ 0.7. The figure shows that, while for low values of x gtegt,FCS is much greater that gte gt,2D, for high values of x the latter becomes greater than the former; however, there is a wide range of values of x, in the neighborhood of x ¼ 1, where the two g-functions are practically coincident. In particular, they coincide for x ¼ 0.8403. This is the optimal value of x for the transition from gtegt,2D to gtegt,FCS, for the case considered. gt,2D

5. Evaluation of the g-functions at the interface tubes-grout, for both low and high values of the dimensionless time, 10L4 £ s* £ 106: results By means of the method illustrated in Section 4, accurate numerical values of 45 g-functions at the interface tubes-grout were obtained, in the range 4  x  6. More precisely, for each value of L*, 9 g-functions were determined, each for a different pair of values of k* and (rc)*. Then, each g-function was plotted versus x ¼ log10s* and interpolated by means of two polynomial functions of x, the first valid for low values of x, the second for high values of x. As in Ref. [34], the interpolating functions were written in the following form:

for  4  x < x0 ; gðxÞ ¼ 0; for x0  x < x1 ; gðxÞ ¼ a6 x6 þ a5 x5 þ a4 x4 þ a3 x3 þ a2 x2 þ a1 x þ a0 ; for x1  x  6; gðxÞ ¼ b6 x6 þ b5 x5 þ b4 x4 þ b3 x3 þ b2 x2 þ b1 x þ b0 :

(13)

Please cite this article in press as: Zanchini E, Lazzari S, New g-functions for the hourly simulation of double U-tube borehole heat exchanger fields, Energy (2014), http://dx.doi.org/10.1016/j.energy.2014.04.022

E. Zanchini, S. Lazzari / Energy xxx (2014) 1e12

7

Table 4 Values of the constants x0, a6, a5, a4, a3, a2, a1, a0, and x1, b6, b5, b4, b3, b2, b1, b0 of the g-function gtegt, for L* ¼ 2000. L* ¼ 2000

(rc)* 0.4 0.7 1 (rc)* 0.4 0.7 1 (rc)* 0.4 0.7 1

x0

a6

a5

a4

a3

a2

a1

a0

x1

b6

b5

b4

b3

b2

b1

b0

0.0034124 e 0.002676 e 0.0022429 e

0.0514849 0.0000975 0.0379012 0.0000928 0.0304264 0.0000788

0.3099394 0.0010183 0.210156 0.0008621 0.1590026 0.0005757

0.9323837 0.0047989 0.5622872 0.0031715 0.3872834 0.0009165

1.3955046 0.0127623 0.6792564 0.0054795 0.3738983 0.0029692

0.762358 0.1610814 0.1148541 0.1753116 0.1173924 0.1902138

0.20705 0.41804 0.37109 0.4098 0.39843 0.40107

0.0013543 e 0.001395 e 0.0013737 e

0.0216449 0.0001228 0.0205703 0.0001223 0.019268 0.0001176

0.1390559 0.0014154 0.1195434 0.0013516 0.1048674 0.001231

0.4475419 0.0070695 0.3366666 0.0061816 0.2682548 0.0050547

0.7083455 0.0184959 0.4256488 0.0139406 0.2762195 0.009054

0.3704648 0.1548869 0.0591066 0.1645771 0.0769356 0.1743766

0.19014 0.32101 0.28027 0.31573 0.30037 0.30935

0.0003251 e 0.0007038 e 0.0008479 e

0.0063209 0.0001279 0.0109134 0.0001296 0.0122353 0.0001263

0.0480842 0.0015093 0.0666421 0.0014774 0.0684911 0.0013854

0.1757904 0.0077191 0.1946881 0.0070081 0.1783495 0.0060905

0.2899713 0.0205875 0.2398516 0.0164653 0.1745852 0.0123143

0.0833979 0.1518138 0.0281073 0.1610414 0.1026249 0.1695997

0.22393 0.28194 0.25784 0.27766 0.26723 0.27216

k* ¼ 0.4 3.996 0.6 3.976 0.35 3.954 0.20 k* ¼ 0.7 4 0.56 3.996 0.3 3.972 0.15 k* ¼ 1 4 0.5 4 0.25 4 0.1

The values of the constants x0, x1, a6, a5, a4, a3, a2, a1, a0, b6, b5, b4, b3, b2, b1, b0, for k* ¼ 1, 0.7, 0.4 and (rc)* ¼ 1, 0.7, 0.4 are reported in Table 4 for L* ¼ 2000, in Table 5 for L* ¼ 1400, in Table 6 for L* ¼ 1000, in Table 7 for L* ¼ 700, in Table 8 for L* ¼ 500. The accuracy of the interpolation is excellent. For instance, for L* ¼ 700, k* ¼ 0.7 and (rc)* ¼ 0.7, which are common values of these dimensionless parameters, the root mean square deviation between the values of gte gt obtained through the interpolating function and those obtained by the numerical computation is 2.13
105; the relative root mean square deviation, i.e., the ratio between the root mean square deviation and the mean value of gtegt, is 4.11
105. When represented in the whole range 4  x < 6, the plots of gtegt(x) obtained by means of Eq. (13) are indistinguishable from those obtained by employing the computational results.

The effect of the thermal conductivity of the grout on the gfunction at the interface tubes-grout, for double U-tube BHEs, is illustrated in Figs. 9 and 10. In Fig. 9, plots of gtegt versus x, in the range 4  x  6, are reported, for L* ¼ 1000 and (rc)* ¼ 0.7, for three different values of k*. The figure shows that low values of the dimensionless thermal conductivity of the grout, such as k* ¼ 0.4, should be avoided because they reduce significantly the thermal performance of the BHEs. The value k* ¼ 0.4 corresponds, for instance, to a cement e bentonite grout with thermal conductivity kgt ¼ 1 W/(m K) and a ground with thermal conductivity kg ¼ 2.5 W/(m K). Similar results, for a different ratio between BHE length and diameter, are illustrated in Fig. 10, where plots of gtegt versus x, in the range 4  x  6, are reported for L* ¼ 700 and (rc)* ¼ 0.7, for different values of k*.

Table 5 Values of the constants x0, a6, a5, a4, a3, a2, a1, a0, and x1, b6, b5, b4, b3, b2, b1, b0 of the g-function gtegt, for L* ¼ 1400. L* ¼ 1400

(rc)* 0.4 0.7 1 (rc)* 0.4 0.7 1 (rc)* 0.4 0.7 1

x0

a6

a5

a4

a3

a2

a1

a0

x1

b6

b5

b4

b3

b2

b1

b0

0.003379 0.0000116 0.002574 0.0000078 0.0021078 0.0000099

0.0509948 0.0000993 0.0364752 0.0000551 0.0285958 0.0001126

0.3070554 0.0002873 0.2022711 0.0002797 0.1492643 0.0009223

0.9237373 0.0011072 0.5404199 0.0007017 0.361601 0.0043454

1.3816703 0.0087607 0.6475799 0.000209 0.3391267 0.0112621

0.7512438 0.1613909 0.0924893 0.1777801 0.1397354 0.194745

0.21053 0.42062 0.37699 0.41141 0.40357 0.4021

0.0016343 0.000001 0.0014793 0.0000069 0.0014021 0.0000104

0.0257283 0.000131 0.0217348 0.0002262 0.0196492 0.0002813

0.1628617 0.0013081 0.1258938 0.0018066 0.1068652 0.0020678

0.5181041 0.0063436 0.3539672 0.0068805 0.2734197 0.006847

0.8196565 0.0170697 0.4501255 0.0138744 0.2830221 0.0103063

0.4583167 0.1554023 0.07584 0.1650122 0.0727371 0.1743845

0.16348 0.32293 0.27606 0.31655 0.29947 0.31001

0.0005026 0.000001 0.0008034 0.0000087 0.0009071 0.0000127

0.0088729 0.0001012 0.0122773 0.0002636 0.0130205 0.0003323

0.0627032 0.001175 0.0739875 0.0021228 0.0725545 0.0025045

0.2181876 0.0063255 0.2143764 0.0082418 0.1886777 0.0087592

0.3550382 0.0183173 0.2671036 0.0170237 0.1878656 0.0148193

0.1329763 0.152642 0.0100362 0.1612901 0.0947191 0.1689242

0.20952 0.28344 0.25329 0.27802 0.26561 0.2727

k* ¼ 0.4 4 0.65 3.974 0.4 3.958 0.25 k* ¼ 0.7 4 0.65 3.99 0.35 3.978 0.2 k* ¼ 1 4 0.55 4 0.3 3.99 0.15

Please cite this article in press as: Zanchini E, Lazzari S, New g-functions for the hourly simulation of double U-tube borehole heat exchanger fields, Energy (2014), http://dx.doi.org/10.1016/j.energy.2014.04.022

8

E. Zanchini, S. Lazzari / Energy xxx (2014) 1e12

Table 6 Values of the constants x0, a6, a5, a4, a3, a2, a1, a0, and x1, b6, b5, b4, b3, b2, b1, b0 of the g-function gtegt, for L* ¼ 1000. L* ¼ 1000

(rc)* 0.4 0.7 1 (rc)* 0.4 0.7 1 (rc)* 0.4 0.7 1

x0

a6

a5

a4

a3

a2

a1

a0

x1

b6

b5

b4

b3

b2

b1

b0

0.003264 0.0000045 0.0024198 0.000012 0.0019387 0.0000142

0.049292 0.000127 0.0343007 0.0002339 0.0262821 0.0002511

0.2969489 0.0006854 0.1901051 0.0010878 0.136798 0.000924

0.8931179 0.0025018 0.5061698 0.0017816 0.3281696 0.0005105

1.3320411 0.0084955 0.596984 0.0006873 0.2928191 0.0088545

0.710722 0.162047 0.0558157 0.1783044 0.1704734 0.1949732

0.22346 0.41953 0.38701 0.41019 0.411 0.40076

0.0017128 0.000018 0.0015176 0.0000303 0.0013925 0.0000386

0.0268894 0.0003726 0.0222698 0.0005799 0.0195188 0.0007203

0.1697441 0.0023852 0.128845 0.003611 0.1061686 0.0044412

0.538922 0.008068 0.3621278 0.0107201 0.2715725 0.0124112

0.8533324 0.0172692 0.4619011 0.0166426 0.2805023 0.0153476

0.4857449 0.1558626 0.0841107 0.1646669 0.0743738 0.1730674

0.15471 0.32225 0.27389 0.31614 0.29985 0.30969

0.000649 0.0000174 0.0008816 0.000033 0.0009454 0.0000421

0.0109979 0.0003663 0.013358 0.0006346 0.0135334 0.0007933

0.07502 0.0023999 0.0798819 0.0040556 0.0752443 0.0050439

0.2544326 0.0084568 0.2304387 0.0125402 0.1956379 0.0149343

0.4116964 0.0189538 0.2898354 0.02058 0.1970442 0.0209665

0.1771778 0.1530369 0.0055113 0.160409 0.0890436 0.1667575

0.19632 0.28327 0.24945 0.27762 0.26439 0.27251

k* ¼ 0.4 4 0.70 3.978 0.45 3.968 0.3 k* ¼ 0.7 4 0.70 4 0.4 3.98 0.25 k* ¼ 1 4 0.6 4 0.35 3.982 0.2

6. Polynomial expressions of the g-functions for L* [ 500 Polynomial expressions of the g-functions at r* ¼ 0.5, 30, 40, 60, 80, 120, 170, 230, 300, 400, 600, for L* ¼ 2000, 1400, 1000, 700, were determined in Ref. [34]; the values of the coefficients of the polynomials are reported there in Tables 4e7. Here, the results obtained in Ref. [34] are extended by providing a table of the coefficients of the polynomial expressions of the g-functions at r* ¼ 0.5, 30, 40, 60, 80, 120, 170, 230, 300, 400, 600, for L* ¼ 500. The g-functions for L* ¼ 500 were determined by the same method of Ref. [34], namely by 2-D axisymmetric numerical simulations performed through COMSOL Multiphysics, where the BHE is sketched as a cylindrical heat source, and polynomial interpolations of the results. A computational domain with dimensionless radius equal to 2500

and dimensionless length equal to 2500 (i.e., L* þ 2000) was employed, together with an unstructured mesh with 33675 triangular elements. The values of the coefficients of the g-functions at r* ¼ 0.5, 30, 40, 60, 80, 120, 170, 230, 300, 400, 600, for L* ¼ 500, are reported in Table 9. 7. Example In this Section, the errors which can be caused by the employment of the quasi-stationary approximation of the heat transfer in the grout, for the hourly simulation of double U-tube BHE fields, are illustrated by an example. Since the internal structure of the BHEs has no relevant effects on the interference between them, only one double U-tube BHE is considered, with L* ¼ 700, placed in a ground

Table 7 Values of the constants x0, a6, a5, a4, a3, a2, a1, a0, and x1, b6, b5, b4, b3, b2, b1, b0 of the g-function gtegt, for L* ¼ 700. L* ¼ 700

(rc)* 0.4 0.7 1 (rc)* 0.4 0.7 1 (rc)* 0.4 0.7 1

x0

a6

a5

a4

a3

a2

a1

a0

x1

b6

b5

b4

b3

b2

b1

b0

0.003264 0.0000346 0.0022123 0.0000413 0.0017345 0.0000457

0.049292 0.000528 0.0313462 0.0006201 0.0234589 0.0006732

0.2969489 0.0024116 0.1733739 0.002703 0.121393 0.0027683

0.8931179 0.0052849 0.4583413 0.0041645 0.2861686 0.0025176

1.3320411 0.0089453 0.5249199 0.0004402 0.2333355 0.0083905

0.710722 0.1630263 0.002203 0.179787 0.2111972 0.1963884

0.22346 0.41857 0.40219 0.4087 0.42128 0.39913

0.003264 0.0000346 0.0022123 0.0000413 0.0017345 0.0000457

0.049292 0.000528 0.0313462 0.0006201 0.0234589 0.0006732

0.2969489 0.0024116 0.1733739 0.002703 0.121393 0.0027683

0.8931179 0.0052849 0.4583413 0.0041645 0.2861686 0.0025176

1.3320411 0.0089453 0.5249199 0.0004402 0.2333355 0.0083905

0.710722 0.1630263 0.002203 0.179787 0.2111972 0.1963884

0.22346 0.41857 0.40219 0.4087 0.42128 0.39913

0.0007617 0.0000471 0.0009321 0.0000671 0.0009584 0.0000788

0.0126485 0.0007601 0.014062 0.0011072 0.0137092 0.0013134

0.0846957 0.0040784 0.0837686 0.0062574 0.0761768 0.0075766

0.2833067 0.011097 0.2411961 0.0167646 0.1980877 0.0202369

0.4576294 0.0192187 0.3053777 0.0229313 0.2003433 0.0249428

0.2138144 0.1540884 0.016446 0.160643 0.0869401 0.1660737

0.18514 0.28303 0.24686 0.27735 0.26391 0.27208

k* ¼ 0.4 4 0.70 3.986 0.5 3.97 0.35 k* ¼ 0.7 4 0.70 3.986 0.5 3.97 0.35 k* ¼ 1 4 0.65 4 0.4 3.982 0.25

Please cite this article in press as: Zanchini E, Lazzari S, New g-functions for the hourly simulation of double U-tube borehole heat exchanger fields, Energy (2014), http://dx.doi.org/10.1016/j.energy.2014.04.022

E. Zanchini, S. Lazzari / Energy xxx (2014) 1e12

9

Table 8 Values of the constants x0, a6, a5, a4, a3, a2, a1, a0, and x1, b6, b5, b4, b3, b2, b1, b0 of the g-function gtegt, for L* ¼ 500. L* ¼ 500

(rc)* 0.4 0.7 1 (rc)* 0.4 0.7 1 (rc)* 0.4 0.7 1

x0

a6

a5

a4

a3

a2

a1

a0

x1

b6

b5

b4

b3

b2

b1

b0

0.003264 0.0000621 0.0020046 0.0000642 0.0015379 0.0000684

0.049292 0.0008565 0.0283594 0.0008652 0.0207135 0.0009131

0.2969489 0.0036062 0.1562623 0.0033277 0.1062257 0.0033256

0.8931179 0.0067018 0.4087086 0.0037757 0.2441523 0.0017588

1.3320411 0.0083141 0.4487358 0.0027718 0.1725616 0.012537

0.710722 0.1635548 þ0.0558589 0.1815621 0.2540204 0.1991266

0.22346 0.41842 0.41916 0.40767 0.43255 0.39808

0.0017128 0.0000756 0.0014758 0.0000885 0.0012854 0.000099

0.0268894 0.0011031 0.0216729 0.0013204 0.0180425 0.0014988

0.1697441 0.0053112 0.1254591 0.0065981 0.0981453 0.0076591

0.538922 0.0122777 0.3524282 0.0150871 0.249812 0.0173707

0.8533324 0.0170872 0.4472473 0.0165819 0.2498979 0.0157955

0.4857449 0.1573824 0.0731712 0.1662939 0.0951276 0.1748108

0.1547 0.32115 0.277 0.31536 0.30502 0.3089

0.0008468 0.0000728 0.0009587 0.0009523 0.0009523 0.0001044

0.0139078 0.0010585 0.0144367 0.0136261 0.0136261 0.0016061

0.0921626 0.0050758 0.085863 0.0757279 0.0757279 0.0085031

0.3059017 0.0119286 0.2470851 0.1968793 0.1968793 0.0206907

0.4941986 0.0178485 0.3140631 0.198661 0.198661 0.0226412

0.2436183 0.1548565 0.0227278 0.0880647 0.0880647 0.1677826

0.17553 0.28311 0.24523 0.26418 0.26418 0.27231

k* ¼ 0.4 4 0.70 3.984 0.55 3.972 0.4 k* ¼ 0.7 4 0.70 3.992 0.5 3.98 0.35 k* ¼ 1 4 0.7 4 0.45 3.99 0.3

with thermal conductivity kg ¼ 2.5 W/(m K), heat capacity per unit volume (rc)g ¼ 2.3 MJ/(m3K) and undisturbed ground temperature Tg ¼ 14  C. A grout with thermal conductivity kgt ¼ 1.61 W/(mK) is considered, so that k* ¼ 0.7. If D ¼ 0.15 m, the BHE length is 105 m. The value of (rc)* is assumed equal to 0.7. The reference heat load per unit length Q0 is considered as equal to the magnitude of the mean monthly heat load per unit length of January, with value Q0 ¼ 30 W/m. The BHE is assumed as working in heating mode during the winter season in the climate of Bologna, North-Center Italy. A stepwise constant heat load is considered, with steps of 1 h; for the given values of kg and (rc)g, the dimensionless duration of 1 h is s*0 ¼ 0.1558. For each hour, the value of the dimensionless heat load per unit length Q* ¼ Q/Q0 is assumed as equal to the ratio between the hourly energy need for heating and the mean value of the magnitude of that quantity during January, for a one family house located in Bologna, determined through TRNSYS 16 [36]. To * avoid too
high
values of the magnitude of Q , the highest allowed value of
Q *
is set equal to 2.5, and the missing energy need is supplied during the following hour. The analytical expression of Q* as a function of the dimensionless time s* is

4391   h    i X Q * s* ¼ Aiþ1 H s*  is*0  H s*  ði þ 1Þs*0 ;

Fig. 9. Plots of gtegt versus x, in the range 4  x  6, with L* ¼ 1000 and (rc)* ¼ 0.7, for k* ¼ 1, 0.7, 0.4.

(14)

i¼0

where 4392 h are considered, from October 15 to April 15, Ai is the value of Q* during the i-th hour (negative for every
value
of i), and H is the Heaviside unit-step function. A plot of
Q *
versus s* =s*0 (number of hours) during the heating season, for the case considered, is reported in Fig. 11. The mean temperature at the interface tubes-grout, Ttegt, is evaluated by two different methods: through the g-function at the interface tubes-grout, gtegt, for L* ¼ 700, k* ¼ 0.7 and (rc)* ¼ 0.7; through the g-function at the interface BHE ground g0.5, for L* ¼ 700, determined in Ref. [34]. The thermal resistance of the grout has been evaluated by a 2-D steady-state simulation performed through COMSOL Multiphysics, and is Rgt ¼ 0.0399 m K/W. The expressions of Ttegt as a function of the dimensionless time s* are, respectively

Fig. 10. Plots of gtegt versus x, in the range 4  x  6, with L* ¼ 700 and (rc)* ¼ 0.7, for k* ¼ 1, 0.7, 0.4.

Please cite this article in press as: Zanchini E, Lazzari S, New g-functions for the hourly simulation of double U-tube borehole heat exchanger fields, Energy (2014), http://dx.doi.org/10.1016/j.energy.2014.04.022

10

E. Zanchini, S. Lazzari / Energy xxx (2014) 1e12

Table 9 Values of the constants x0, a6, a5, a4, a3, a2, a1, a0, and x1, b6, b5, b4, b3, b2, b1, b0 of the g-functions at r* ¼ 0.5, 30, 40, 60, 80, 120, 170, 230, 300, 400, 600, for L* ¼ 500. L* ¼ 500 r* 0.5 30 40 60 80 120 170 230 300 400 600

x0 3.826 1.458 1.748 2.184 2.414 2.672 3.116 3.372 3.61 3.79 4.102

a6 0.000068 0.014241 0.0064649 0.0072125 0.0103164 0.0146456 e e e 0.01092804 e

a5 0.0004271 0.195288 0.0997785 0.1248228 0.1961421 0.3033699 e 0.0101022 0.0063342 0.2904565 0.00360603

a4 0.0007473 1.079672 0.6124448 0.8692397 1.5189119 2.5802257 0.0302216 0.2273364 0.1563745 3.18835846 0.07326077

a3 0.0061278 3.062256 1.8967375 3.1008077 6.1234437 11.5255014 0.4490819 2.0084646 1.4988034 18.508177 0.58416854

a2

r* 0.5 30 40 60 80 120 170 230 300 400 600

x1 1.35 3.15 3.6 3.93 4.21 4.3 4.4 4.77 4.91 4.8 5.13

b6 0.0001913 0.0022582 0.0039354 e e e e e e e e

b5 0.0034131 0.0625441 0.1105962 0.0069936 0.004601 0.0060394 0.0085694 0.0234604 e 0.02101006 0.01326315

b4 0.0233535 0.7084645 1.2782469 0.1801604 0.1196875 0.1583825 0.2263958 0.6284542 0.0100843 0.56733137 0.37404103

b3 0.0818113 4.1978044 7.775852 1.8344654 1.2262057 1.6424229 2.373444 6.71044 0.2321316 6.1109557 4.21824322

b2

 

Ttgt s*

h   X Q 4391 ¼ Tg þ 0 Aiþ1 gtgt s*  is*0 kg i¼0  i  gtgt s*  ði þ 1Þs*0 ;

(15)

  h    i X Q 4391 Ttgt s* ¼ Tg þ 0 Aiþ1 g0:5 s*  is*0  g0:5 s*  ði þ 1Þs*0 kg i ¼ 0   þ Q0 Rgt Q * s* : (16)

s* =s*0

Two plots of Ttegt versus for the whole heating season, obtained through Eqs. (15)e(16) respectively, are reported in Fig. 12. The plot obtained through Eq. (15), i.e. by employing the new g-function at the interface tubes-grout, is drawn in black; that obtained through Eq. (16) is drawn in orange (light gray when printed in black and white). A plot of the difference between the value of Ttegt obtained through Eq. (16) and that obtained through Eq. (15) is reported in Fig. 13. The figures show that the discrepancy



Fig. 11. Plot of
Q *
versus s* =s*0 for the example considered.

0.0221063 4.715197 3.1285719 5.9715622 13.5714625 28.5308593 2.4544879 8.7053136 6.9977183 59.95404 2.27894721 0.1538121 13.7461932 26.291753 9.1972591 6.1549494 8.3928008 12.3203117 35.683985 2.0066071 32.8098137 23.7816547

a1

a0

0.1477453 3.755953 2.6082876 5.8818568 15.709289 37.1454893 5.8758636 18.5362496 15.976323 102.815066 4.32926784

0.1955375 1.214788 0.8565059 2.3078553 7.4338477 19.8916775 5.2151692 15.5406194 14.317663 72.9655685 3.1784414

b1

b0

0.0308441 23.7790033 47.0563098 22.584824 15.0212218 21.0365441 31.5859797 94.4436902 7.723773 87.7648196 67.0375427

0.23595 17.04851 34.93571 21.75163 14.2404 20.6612 31.960534 99.504275 11.134982 93.540235 75.595094

between the two plots is neither very high nor quite negligible, and that the quasi-stationary approximation yields an overestimation of Ttegt with this kind of heat load. Plots of Ttegt obtained through Eq. (16) and through Eq. (15) in the narrower range 2180  s* =s*0  2250 (70 h during January) are reported in Fig. 14, and a plot of the difference between the two functions, in the same time interval, is reported in Fig. 15. The mean value of the difference, in this interval, is 0.653  C; much higher values, which can exceed 2.0  C, occur during a few hours. Finally, to verify that the differences illustrated in Figs. 12e15, not predictable intuitively, are not due to inaccuracies of the gfunctions gtegt and g0.5, the comparison between the values of Ttegt obtained through Eq. (15) and those obtained through Eq. (16) was repeated by using a more regular dimensionless heat load, with value 1 during odd hours and value þ1 during even hours. The comparison is illustrated in Fig. 16. In this test case, the results are as expected: the quasi-stationary approximation for the heat transfer in the grout yields a lower value of Ttegt when the heat load

Fig. 12. Plots of Ttegt versus s* =s*0 for the whole heating season, obtained through Eq. (15) (black line) and through Eq. (16) (orange line). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Please cite this article in press as: Zanchini E, Lazzari S, New g-functions for the hourly simulation of double U-tube borehole heat exchanger fields, Energy (2014), http://dx.doi.org/10.1016/j.energy.2014.04.022

E. Zanchini, S. Lazzari / Energy xxx (2014) 1e12

Fig. 13. Difference between the value of Ttegt obtained through Eq. (16) and that obtained through Eq. (15), for the whole heating season.

11

Fig. 16. Difference between the value of Ttegt obtained through Eq. (16) and that obtained through Eq. (15), if Q* ¼ 1 during odd hours and Q* ¼ 1 during even hours, for 0  s* =s*0  20.

recommended for the simulation of double U-tube BHEs. However, the errors are not very high, at least in the case examined; therefore, the use of accurate g-functions at the interface BHE-ground, combined with the quasi-stationary approximation within the BHE, can be considered as acceptable for different BHE geometries, in the absence of g-functions at more internal surfaces. 8. Conclusions

Fig. 14. Plots of Ttegt versus s* =s*0 obtained through Eq. (15) (black line) and through Eq. (16) (orange line), in the range 2180  s* =s*0  2250. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

is negative and a higher value when the heat load is positive; moreover there is no accumulation of the errors when time goes by. The results obtained allow to conclude that the quasi-stationary approximation for the heat transfer in the grout can cause nonnegligible errors in the hourly simulation of BHE fields, so that the use of the new g-functions at the interface tubes-grout is

Polynomial expressions of new g-functions have been determined, which yield the dimensionless temperature at the interface tubes-grout for double U-tube BHEs; the new g-functions allow more accurate hourly simulations of this kind of BHEs with respect to the usual g-functions obtained by the finite-cylindrical-source model or by the finite-line-source model. Tables of the coefficients of the new g-functions have been reported, for several values of the ratio between the BHE length L and the BHE diameter D, and of the ratios between the thermal properties of the ground and those of the grout. With reference to a specific case, the results obtained by the new g-functions at the interface tubes-grout have been compared with those obtained through the g-functions at the interface BHE-ground determined in (Zanchini E., Lazzari S., Energy 2013; 59: 570e80). The results show that the quasi-stationary approximation for the heat transfer in the grout can cause non-negligible errors, so that the use of the new g-functions is recommended for the hourly simulation of double U-tube BHEs. Acknowledgment This research was funded by the European Union’s Seventh Framework Program, Theme EeB.NMP.2012-2, Project HERB (Holistic energy-efficient retrofitting of residential buildings), Grant agreement no: 314283. References

Fig. 15. Difference between the value of Ttegt obtained through Eq. (16) and that obtained through Eq. (15), in the range 2180  s* =s*0  2250.

[1] Chap. 32 ASHRAE handbook e HVAC applications; 2007. [2] Gehlin S. Thermal response test: method development and evaluation. Doctoral thesis. Luleå: Luleå University of Technology; 2002. [3] Rainieri S, Bozzoli F, Pagliarini G. Modeling approaches applied to the thermal response test: a critical review of the literature. HVAC&R Res 2011;17:977e 90. [4] Bozzoli F, Pagliarini G, Rainieri S, Schiavi L. Estimation of soil and grout thermal properties through a TSPEP (two-step parameter estimation procedure) applied to TRT (thermal response test). Energy 2011;36:839e46. [5] Gehlin S, Hellström G. Comparison of four models for thermal response test evaluation. ASHRAE Trans 2003;109:131e42.

Please cite this article in press as: Zanchini E, Lazzari S, New g-functions for the hourly simulation of double U-tube borehole heat exchanger fields, Energy (2014), http://dx.doi.org/10.1016/j.energy.2014.04.022

12

E. Zanchini, S. Lazzari / Energy xxx (2014) 1e12

[6] Wagner R, Clauser C. Evaluating thermal response tests using parameter estimation for thermal conductivity and thermal capacity. J Geophys Eng 2005;2:349e56. [7] Bujok P, et al. Assessment of the influence of shortening the duration of TRT (thermal response test) on the precision of measured values. Energy 2014;64: 120e9. [8] Yavuzturk C, Spitler JD, Rees SJ. A transient two-dimensional finite volume model for the simulation of vertical U-tube ground heat exchangers. ASHRAE Trans 1999;105:465e74. [9] Zeng H, Diao N, Fang Z. Heat transfer analysis of boreholes in vertical ground heat exchangers. Int J Heat Mass Transf 2003;46:4467e81. [10] Lamarche L, Beauchamp B. New solutions for the short-time analysis of geothermal vertical boreholes. Int J Heat Mass Transf 2007;50:1408e19. [11] Zanchini E, Lazzari S, Priarone A. Effects of flow direction and thermal shortcircuiting on the performance of small coaxial ground heat exchangers. Renew Energy 2010;35:1255e65. [12] Zanchini E, Lazzari S, Priarone A. Improving the thermal performance of coaxial borehole heat exchangers. Energy 2010;35:657e66. [13] Rybach L, Eugster WJ. Sustainability aspects of geothermal heat pump operation, with experience from Switzerland. Geothermics 2010;39:365e9. [14] Lazzari S, Priarone A, Zanchini E. Long-term performance of BHE (borehole heat exchanger) fields with negligible groundwater movement. Energy 2010;35:4966e74. [15] Chiasson AD, Rees SJ, Spitler JD. A preliminary assessment of the effects of groundwater flow on closed-loop ground-source heat pump systems. ASHRAE Trans 2000;106:380e93. [16] Zanchini E, Lazzari S, Priarone A. Long-term performance of large borehole heat exchanger fields with unbalanced seasonal loads and groundwater flow. Energy 2012;38:66e77. [17] Lee CK, Lam HN. A modified multi-ground-layer model for borehole ground heat exchangers with an inhomogeneous groundwater flow. Energy 2012;47: 378e87. [18] Li M, Lai ACK. New temperature response functions (G functions) for pile and borehole ground heat exchangers based on composite-medium line-source theory. Energy 2012;38:255e63. [19] Loveridge F, Powrie W. Temperature response functions (G-functions) for single pile heat exchangers. Energy 2013;57:554e64. [20] Loveridge F, Powrie W. G-Functions for multiple interacting pile heat exchangers. Energy 2014;64:747e57.

[21] Eskilson P. Thermal analysis of heat extraction boreholes. PhD thesis. Sweden: Lund University; 1987. [22] Claesson J, Eskilson P. Conductive heat extraction by a deep borehole, analytical studies. Tech. Rep. Sweden: Lund University; 1987 [23] Kavanaugh S, Rafferty K. Ground-source heat pumps: design of geothermal systems for commercial and industrial buildings. Atlanta, GA: ASHRAE; 1997. [24] Hellström G, Sanner B. EED e earth energy designer, user manual. Version 2.0, http://buildingphysics.com/manuals/eed.pdf; 2000. [25] Zeng H, Diao N, Fang Z. A finite line-source model for boreholes in geothermal heat exchangers. Heat Transf Asian Res 2002;31:558e67. [26] Lamarche L, Beauchamp B. A new contribution to the finite line-source model for geothermal boreholes. Energy Build 2007;39:188e98. [27] Bandos TV, Montero A, Fernandez E, Santander JLG, Isidro JM, Perez J, et al. Finite line-source model for borehole heat exchangers: effect of vertical temperature variations. Geothermics 2009;38:263e70. [28] Fossa M. The temperature penalty approach to the design of borehole heat exchangers for heat pump applications. Energy Build 2011;43:1473e9. [29] Fossa M. A fast method for evaluating the performance of complex arrangements of borehole heat exchangers. HVAC&R Res 2011;17(6):948e58. [30] Bernier MA, Chahla A, Pinel P. Long-term ground-temperature changes in geoexchange systems. ASHRAE Trans 2008;114(2):342e50. [31] Cimmino M, Bernier M, François A. A contribution towards the determination of g-functions using the finite line source. Appl Therm Eng 2013;51: 401e13. [32] Bernier M, Pinel P, Labib R, Paillot R. A multiple load aggregation algorithm for annual hourly simulations of GCHP systems. HVAC&R Res 2004;4:471e87. [33] Fossa M, Minchio F. The effect of borefield geometry and ground thermal load profile on hourly thermal response of geothermal heat pump systems. Energy 2013;51:323e9. [34] Zanchini E, Lazzari S. Temperature distribution in a field of long Borehole Heat Exchangers (BHEs) subjected to a monthly averaged heat flux. Energy 2013;59:570e80. [35] Carlslaw HS, Jaeger JC. Conduction of heat in solids. Oxford: Oxford University Press; 1959. [36] Zanchini E, Lazzari S. Design of a nearly zero energy one-family house in North-Centre Italy. Vancouver, Canada. In: 11th International Conference on Sustainable Energy Technologies (SET 2012); 2e5 September 2012. pp. 1609e21.

Please cite this article in press as: Zanchini E, Lazzari S, New g-functions for the hourly simulation of double U-tube borehole heat exchanger fields, Energy (2014), http://dx.doi.org/10.1016/j.energy.2014.04.022