Numerical evaluation of the effects of groundwater flow on borehole heat exchanger arrays

Renewable Energy 52 (2013) 230e240

Contents lists available at SciVerse ScienceDirect

Renewable Energy journal homepage: www.elsevier.com/locate/renene

Numerical evaluation of the effects of groundwater flow on borehole heat exchanger arrays Jung Chan Choi a, *, Joonsang Park a, Seung Rae Lee b a b

Norwegian Geotechnical Institute, P.O. Box 3930 Ullevål Stadion, N-0806 Oslo, Norway Department of Civil and Environmental Engineering, KAIST, Daejeon 305-701, Republic of Korea

a r t i c l e i n f o

a b s t r a c t

Article history: Received 6 June 2012 Accepted 17 October 2012 Available online 24 November 2012

The effects of both direction and rate of groundwater flow on the performance of various types of borehole heat exchanger (BHE) arrays are examined using a two-dimensional coupled heat conduction eadvection model. The heating operations were simulated over a period of 15 yrs using three types of arrays: rectangular, L-type, and single line. The results show that the performance of the L-type and single line-type arrays (non-square rectangular arrays) was noticeably influenced by the direction of groundwater flow as well as the flow rate. When the characteristic length of Peclet number was assumed to be a unit value, the Peclet number less than 0.05 (i.e. low flow rates less than 1 m/yr in this study) was found to have negligible effects on the performance of the system, regardless of the array type or flow direction. The cold accumulation in the ground downstream of the groundwater flow seems to be related to the variation of the fluid temperature associated with the flow direction. The comparison of annual heat capacity shows that up to 13% difference can occur depending on the flow direction. This suggests that the consideration of both direction and rate of groundwater flow may be important in designing the optimal BHE arrays, particularly for the non-square rectangular arrays. Ó 2012 Elsevier Ltd. All rights reserved.

Keywords: Geothermal Borehole heat exchanger array Groundwater flow Coupled heat conductioneadvection model

1. Introduction A ground source heat pump (GSHP) system has received much attention as a renewable energy resource for space heat and cooling due to its high energy efficiency and environment-friendly benefits [1e4]. A vertical GSHP system generally consists of a heat pump and a set of multiple borehole heat exchangers (BHEs), which are also called vertical ground heat exchangers (GHEs), as shown in Fig. 1. The heat pump is a mechanical device that transfers heat energy from the ground to a building, and vice versa; its performance mainly depends on the temperature of the circulating fluid in the BHEs. Hence, the accurate evaluation of heat transfer between BHEs and the surrounding ground is important for an optimal design of a ground source heat pump system. To estimate the optimal size of BHEs for the direct use of geothermal energy, many approaches have been suggested via analytical or numerical methods. The infinite line source model [5] and infinite cylindrical source model [6] are the simplest ones based on transient analytical solutions. These analytical solutions regard the ground as an infinite medium and simplify heat transfer

* Corresponding author. Tel.: þ47 90 29 16 61; fax: þ47 22 23 04 48. E-mail address: [email protected] (J.C. Choi). 0960-1481/$ e see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.renene.2012.10.028

as a one-dimensional conduction in the perpendicular direction to the borehole axis. Eskilson [7] proposed a semi-analytical solution of the finite line source model to account for the finite length and configuration of the borehole using the thermal response function (so called g-function). Zeng et al. [8] and Lamarche and Beauchamp [9] also presented an analytical solution for the finite line source model. Hellström [10] developed a duct storage model (DST) to calculate the performance of a thermal energy storage, which assumes a closely packed layout of boreholes in a cylindrical storage. Yavuzturk and Spitler [11] examined the transient short time response inside a borehole using the two-dimensional finite volume method. Recently, to simulate the complex thermal behavior of BHEs and their vicinity more accurately, numerical approaches based on the finite element (FE) method have also been suggested and the results were compared with experimental data [12e14]. Most previous approaches considered only conductive heat flow around the boreholes and they provide a good estimation of BHE systems. However, when there is groundwater flow around the boreholes, heat flux between BHEs and the surrounding ground may be underestimated due to an ignorance of groundwater advection [15]. Heat advection by groundwater flow is induced by a hydraulic gradient in the ground and highly influenced by the permeability. The advection can affect the performance of a BHE system

J.C. Choi et al. / Renewable Energy 52 (2013) 230e240

Nomenclature

Symbol a c g h hb H K k n Pr pw Re r R T

u

thermal diffusivity (m2 s1) heat capacity (J kg1 K1) acceleration of gravity (m s2) convective heat transfer coefficient of the circulating fluid (W m2 K1) equivalent convective heat transfer coefficient inside the borehole (W m2 K1) depth from groundwater level (m) hydraulic conductivity (m s1) intrinsic permeability (m2) unit normal vector () Prandtl number () pore pressure (Pa) Reynolds number of the circulating fluid () radius (m) thermal resistance (m K W1) temperature (K or ◦C)

positively and few studies investigating this topic have been carried out recently. Chiasson et al. [16] reported the effects of groundwater flow on the thermal performance of single and 4  4 boreholes in various geological conditions using a numerical approach. An analytical solution which can consider groundwater advection as a conduction problem with a moving infinite line heat source has been introduced [6,17,18]. Fan et al. [19] examined the effects of advection on the underground heat/cold storage by a numerical approach. Also, some numerical investigations on the influence of groundwater flow on the thermal/geothermal response test (TRT/ GRT) have been performed [20e22]. However, most of the researches considering effects of advection only focus on the rate of groundwater flow, but not on the direction. In the case of noncircular or non-rectangular arrays, the direction of groundwater flow may also affect heat transfer in the vicinity of BHEs. Some approaches to evaluate the effects of flow direction on the performance of group BHE systems have been performed [23]. However,

l m r ɸ

231

groundwater flow velocity vector (m s1) thermal conductivity (W m1 K1) viscosity of groundwater (Pa s) density (kg m3) soil porosity ()

Subscripts b borehole eq equivalent pipe f circulating fluid g ground grout grout/filling material in inlet of pipe out outlet of pipe p U-shaped pipe s soil particles w groundwater Superscripts i inner parts of pipe/equivalent pipe o outer parts of pipe/equivalent pipe

the investigation on the various type of non-circular or -rectangular types of BHE arrays is lacking. In this study, the effects of groundwater flow on BHE arrays were evaluated using a coupled heat conductioneadvection transient model. Heat transfer inside the borehole was treated as steady-state thermal resistance and heat transfers outside the boreholes were simulated using a 2-D transient finite element method. The suggested numerical model was compared with a finite line source model for validation. The effects of direction and rate of groundwater flow on the long-term performance of three different types of BHE arrays were examined. 2. Computational model The thermal behavior of a BHE array is complicated because of the various heat transfers inside and outside the BHEs. Heat transfer inside the borehole can be decomposed into the following

Fig. 1. A schematic diagram of the ground source heat pump (GSHP) system.

232

J.C. Choi et al. / Renewable Energy 52 (2013) 230e240

mechanisms: forced convection of circulating fluid inside the Upipe; conductive heat transfer in the pipe wall and the borehole filling material; and thermally induced natural convection (in the case of a water filled borehole). Heat transfer outside the borehole, i.e. ground heat transfer, can be mainly described with heat conduction and advection if there is no groundwater freezing or evaporating. Full implementation of the aforementioned various heat transfer mechanisms into a numerical simulation may explain thermal behavior of a ground heat exchanger system in detail, however a simulation with full implementation of each heat transfer mechanism may be time consuming and costly, particularly for long-term simulations or large borehole arrays [24,25]. As the main focus for this study was the investigation of long-term behavior of BHE arrays affected by groundwater flow rather than a detailed examination of inside-borehole heat transfer, the heat transfer inside the borehole was assumed to be steady-state thermal resistance, and heat transfer outside the borehole was simulated by 2-D transient numerical approaches. As the model simulated the problem in 2-D plane, a constant heat flux along the vertical dimension was assumed. In addition, peak heating demands was not considered in this study. Heat transfer in the ground occurs mainly by conduction in solid soil particles and liquid water phase, and advection through moving liquid water. Although latent heat transfer caused by phase change (e.g., freezing or evaporation) could occur in the ground, it was assumed to be negligible in this study for simplicity. If the local temperatures of groundwater and soil particles, which can be expressed as Tw and Ts respectively, are assumed to reach instantly in thermal equilibrium in fully saturated ground, the governing equation for the coupled heat conduction and advection can be expressed as [19,26]:

rg cg

dT þ rw cw uVT  VðlVTÞ ¼ 0 dt

(1)

where r is the mass density, c is the heat capacity, l is the thermal conductivity of the ground, u is the groundwater flow velocity vector, and T (¼Tw ¼ Ts) is the ground temperature. The subscripts g and w indicate the ground and the groundwater, respectively. The volumetric heat capacity of the ground can be presented as the parallel sum of the thermal properties of each phase according to their volume fraction. However, the thermal conductivity is more influenced by the complex microstructures and able to be estimated by various models [27,28]. The volumetric heat capacity using parallel sum model and the thermal conductivity using weighted geometric mean model [27] can be expressed as:

rg cg ¼ ð1  4Þrs cs þ 4rw cw lg ¼ lð14Þ l4w s

(2)

At the interface between the ground and the borehole, heat exchange can be expressed as following mixed Neumann’s boundary condition:

l


  vT

¼ hb Tf  Tb
vn on borehole wall

where Tf is the mean temperature of the circulating fluid (defined as Tf ¼ (Tf,in þ Tf,out)/2) with subscripts in and out describing the inlet and outlet of pipe; Tb is the borehole wall temperature, which can be assumed to be the same as the ground temperature on the borehole wall; and n is the unit vector normal to the borehole axis. Although Tf calculated by the arithmetic average of Tf,in and Tf,out can be different from the real mean fluid temperature [29], for simplicity, the arithmetic average form was used as a representative of the fluid temperature in this study. hb is the equivalent convective heat transfer coefficient inside the borehole and can be expressed with borehole thermal resistance, Rb, as follows:

hb ¼

1 1 2prb Rb

(6)

where rb is the borehole radius. The borehole thermal resistance, Rb, can generally be obtained by thermal response tests (TRT) or analytical approaches. Although many analytical models for the borehole thermal resistance have been suggested [8,30e32], for simplicity, the one-dimensional series model [30] was adopted in this study. If the U-shaped pipes inside the borehole are approximated to be a single pipe with an equivalent diameter, the borehole thermal resistance can be simplified to be a connection in series of thermal resistances of circulating fluid, Rf, pipe wall, Rp, and grout/ filling material, Rgrout, as shown in Fig. 2, and can be expressed as follows:

Rb ¼ Rf þ Rp þ Rgrout Rf ¼ Rp ¼

Rgrout

1 2prpi h 1 2plp

ln 1

o req

!

i req

r ¼ ln ob req 2plgrout

(7) !

where rp is the U-shaped pffiffiffiffiffiffiffiffi pipe radius, and req is the equivalent pipe radius denoted as 2rp [30]. Superscripts i and o denote the inner and outer parts of pipe/equivalent pipe, respectively. lp is the thermal conductivity of the U-shaped pipes, and lgrout is the thermal conductivity of the grout/filling material. h is the

where ɸ is the soil porosity. The subscripts s indicates the soil particle. The groundwater velocity vector u can be expressed in terms of Darcy’s law as follows:

k u ¼  Vðpw þ rw gHÞ

m

(3)

where k is the intrinsic permeability, m is the viscosity of groundwater, pw is the pore pressure, g is the acceleration of gravity, and H is the depth from groundwater level. If the ground is fully saturated with groundwater, Eq. (3) can also be expressed in terms of hydraulic conductivity, K, as follows:

K ¼

krw g

m

(4)

(5)

Fig. 2. The thermal resistance network inside the borehole.

J.C. Choi et al. / Renewable Energy 52 (2013) 230e240

convective heat transfer coefficient of the circulating fluid and can be calculated from the DittuseBoelter correlation for turbulent flow as follows [33]:

h ¼ 0:023Re0:8 Pr 0:35

lf

(8)

2rpi

where Re is the Reynolds number of the circulating fluid; Pr is the Prandtl number; and lf is the thermal conductivity of the circulating fluid. The above model was integrated into a commercial multiphysics finite element (FE) code COMSOL Multiphysics 4.2a with a convection-diffusion mathematics module. The model was assumed to be a two-dimensional transient heat transfer problem. The ground is assumed to be a homogeneous porous medium fully saturated with groundwater. The effect of the finite length of the borehole (i.e. heat/cold accumulation at the tip of the borehole) was not considered in this study. To avoid spurious oscillations of the calculated results, which generally occur in convectionediffusion type FE calculations, the streamline upwind PetroveGalerkin (SUPG) formulation was adopted for stabilization of the calculation [34]. The domain near the borehole, which shows a steep temperature gradient, was modeled using a fine mesh for accurate calculations. Far field from the pipe, a coarse mesh was generated to reduce the size of memory consumption and the calculation time. 3. Model validation The model was compared with analytical solutions of the moving infinite line source (MILS) model [17,18,32] and Eskilson’s finite line source model (FLSM) with superposition [7] for its validation. The MILS model was used for the validation of the numerical results of the transient heat conductioneadvection analysis and the FLSM with superposition was used for the validation of the numerical results of long-term calculation on the group BHE arrays. For the validation of two-dimensional heat conductione advection numerical model, the calculated temperature response of a single borehole with 50 h simulation was compared with the solution of the MILS model. If the groundwater flows in direction of 4 ¼ 0 with uniform velocity of uw, the analytical solution for the MILS model around the borehole in the polar coordinate can be derived as follows [18]:

 4aZg t=r Ur 1 DTðr; fÞ ¼ exp cos4 exp h 2ag 4plg q



2

0

! U 2 r2 h dh   h 16a2g 1

(9) where q is the heat exchange rate per unit length of the borehole, U ¼ uwrwcw/rgcg is the effective heat transport velocity, and ag is the thermal diffusivity of the ground. Fig. 5 shows the comparison of analytical and numerical solutions. Thermal properties of the ground were selected as rg ¼ 2190 kg/m3, cg ¼ 1050 J/kg/K, and lg ¼ 3.0 W/m/K. The heat exchange rate per unit length and the groundwater flow velocity were selected as q ¼ 10 W/m, and u ¼ 1e5 m/s, respectively. The isotherm contour of temperature change DT at simulation time of 50 h (Fig. 5a) shows good agreement between numerical and analytical solutions. The numerical solution of temperature changes against time (Fig. 5b) also gives good agreement with the analytical solution. The model was also compared with the solution of FLSM with superposition for the validation of its numerical results for the long-term calculation on the group BHE arrays. Eskilson’s model

233

Table 1 Borehole configurations for simulation.

Configuration Borehole spacing (m)

Case 1

Case 2

Case 3

Single line 10.0

L-type 10.0

Rectangular 10.0

Layout

uses a combination of the numerical finite difference method and an analytical solution based on the finite line source model using a temperature response function (known as the g-function). It can consider thermal interferences between boreholes. It has thus been implemented into many commercial ground energy simulation software for calculating the thermal response of BHE arrays although it cannot consider the effects of the groundwater flow. In this study, EED 3.16 [35] was adopted for the finite line source model with superposition. No groundwater flow was considered in the comparison problem. Three different kinds of BHE arrays, with rectangular and nonrectangular shapes, were selected as shown in Table 1. Each array had 9 boreholes to the depth of 200 m; each borehole consisted of a single U-shaped pipe and borehole filling materials. The geological ground condition was assumed to be sedimentary limestone and the thermal properties of the ground were selected to be a value within the representative range of sedimentary limestone [16]. The thermal resistance of the borehole for both numerical and analytical model was calculated using the same method expressed in Eqs. (6)e(8). Additional details about the borehole and the ground for the simulation are summarized in Table 2. In central and northern Europe, the geothermal heat exchanger system is mainly used for heating rather than cooling [37]. The simulation was thus carried out over a period of 15 years, with heating as the prevailing operation, in order to consider long-term effects. Total annual load was assumed to be 180 MWh yr1, considering only heating, which was estimated to be representative of the heat load for small schools or buildings. A seasonal performance factor (SPF) for heating was selected to be a value of 3.0. Monthly distributions of total and ground energy loads are summarized in Fig. 3. The ground

Table 2 Input properties related to the borehole. Parameters Borehole properties Borehole length, L Borehole radius, rb Pipe outer radius, rpo Pipe thickness, tp Distance between centers of pipes, Lshank Pipe thermal conductivity, lp Grout/Filling thermal conductivity, lgrout Effective borehole thermal resistance, Rb Ground properties Ground thermal conductivity, lg Ground density, rg Ground heat capacity, cg Groundwater properties Groundwater thermal conductivity, lw Groundwater density, rw Groundwater heat capacity, cw Circulating fluid properties Fluid thermal conductivity, lf Fluid density, rf Fluid heat capacity, cf Fluid volumetric flow rate, V_

Value

Unit

200.0 0.07 25.0 2.9 70.0 0.42 1.0 0.17

m m mm mm mm W m1 K1 W m1 K1 m K W1

3.0 2190.5 1050.0

W m1 K1 kg m3 J kg1 K1

0.6 1000.0 4180.0

W m1 K1 kg m3 J kg1 K1

0.453 1068.0 3565.0 0.0007

W m1 K1 kg m3 J kg1 K1 m3 s1

234

J.C. Choi et al. / Renewable Energy 52 (2013) 230e240

the finite line source model (FLSM) are shown in Fig. 6 and Table 3, respectively. The results show relatively similar trends for the fluid temperature variation during simulation. However, the numerical solutions resulted in slightly lower mean fluid temperature in longterm simulation, and the rectangular array (Case 3) gave a maximum temperature difference of 0.3  C between the FE model and the FLSM model. This difference appears to be related to the effects of the finite length of boreholes. Consideration of the finite length of boreholes is known to give slightly higher fluid temperature in long-term heat prevailing simulations [38e40]. However, in the scope of this study (i.e. examination of the influence of groundwater flow on BHE arrays), this numerical model gave good agreement with the finite line source model within a relatively small error.

30.0 Total heating load Ground heating load

Energy load (MWh)

25.0 20.0 15.0 10.0 5.0 0.0 JAN FEB MAR APR MAY JUN JUL AUG SEP OCT NOV DEC

4. Results and discussion

Month As briefly mentioned in the Introduction, groundwater flow is known to affect the performance of geothermal heat pump systems positively. It has been suggested that faster groundwater speed increases efficiency of BHE system operation [15,19]. In this section, more detailed examination of various types of BHE arrays is presented; the effects of both the groundwater flow rate and direction on the previously selected three types of BHE arrays were simulated using a coupled conductioneadvection model. The hydraulic conductivity in Eq. (4), which is related to groundwater flow in the coupled heat conduction and advection equation (Eq. (1)), may depend on several factors (e.g., porosity, degree of fracture, mineral composition, etc.). According to the previous research [41,42], the hydraulic conductivity of sedimentary porous limestone roughly ranges within 1.0e7e1.0e4 m/s. In this study, groundwater flow rates ranging from 0 (no flow) to 10 m/yr were thus applied to the simulation, with a constant hydraulic gradient assumption of 0.01 m/m. Steady-state Darcy flow was assumed. The input data about the borehole and the ground were assumed to be the same as in the previous validation problem (which are summarized in Tables 1 and 2). The same energy load in Fig. 3 was also applied to the simulations. The effects of groundwater flow on the maximum drop in mean fluid temperature were examined, as shown in Fig. 6. The maximum drop in mean fluid temperature can be defined by the difference in initial fluid temperature and minimum mean fluid temperature during simulations. Although some energy storage systems propose large decreases in circulating fluid temperature in order to increase storage capacity, a smaller fluid temperature drop

Fig. 3. Monthly distributions of the total and ground energy load.

Fig. 4. The 2-D finite element mesh for Case 2 (L-type configuration).

energy loads were applied to the heat boundary on the boreholes as expressed on the left side in Eq. (5). The initial ground and far-field temperatures were set to 9.0  C. The ground domain was modeled with dimensions of 1000 m (L)  1000 m (W) to consider a suitable margin for the heat front from the heat source. The FE mesh for simulations was comprised of 7718 unstructured mesh with triangular elements, as shown in Fig. 4. The monthly drop and the annual maximum drop in mean fluid temperature calculated by the FE numerical model (this study) and

a

b

2

Analytical solution … … Numerical solution 0.3

Temperature change [K]

1 0.05

0.2 0 15 0. 0.3 .25 5 2 . 0.2 0 0.15 0.1

0.1

0.15

0.1

0.05

0. 05

0.05

y [m]

0.1

05 0.

0.5

-0.5

◦

r=0.5m, φ=0

1.5

0

0.4

-1

Analytical solution Numerical solution 0.2

0.1 ◦

r=1.0m, φ=45

-1.5

t = 50 hrs -2

0.0 0

1

2

x [m]

3

4

0.0

10.0

20.0

30.0

40.0

50.0

Time (hr)

Fig. 5. Comparison of analytical (moving infinite line source model) and numerical solution: a) isotherms contours of temperature change at t ¼ 50 h, and b) temperature changes with time at r ¼ 0.5 m with 4 ¼ 0 and r ¼ 1.0 m with 4 ¼ 45 .

J.C. Choi et al. / Renewable Energy 52 (2013) 230e240

235

Table 3 Annual maximum drop in mean fluid temperature obtained by FLSM and the FE model. Time (yr)

1 2 5 10 15

Fig. 6. Temperature drop in circulating fluid calculated by the finite element (FE) model and the finite line source model (FLSM) for a) single line-type, b) L-type, and c) rectangular borehole arrays.

means more sustainable performance of general heat pump systems, as most heat pumps operate most efficiently within a temperature range of 0e30  C. Excessively low fluid temperatures can also increase uncertainties related with ground freezing around boreholes (i.e. positive effects of latent heat delivers or high thermal conductivity of the freezing ground; negative effects of no heat advection around the borehole or structural instability of the borehole). When there was no groundwater flow simulated (black

Annual maximum drop in mean fluid temperature ( C) Case 1 (Single line-type)

Case 2 (L-type)

Case 3 (Rectangular)

FLSM

FE model

FLSM

FE model

FLSM

FE model

5.4 5.9 6.7 7.3 7.7

5.4 5.9 6.6 7.3 7.8

5.5 5.9 6.7 7.4 7.8

5.5 5.9 6.7 7.5 8.0

5.6 6.1 7.3 8.2 8.8

5.5 6.1 7.3 8.4 9.1

solid rectangular dots in Fig. 7), the mean fluid temperature after a 15-yr simulation dropped to 7.8, 8.0, and 9.1  C for the linetype, L-type, and rectangular arrays, respectively. In most cases, when groundwater flow was considered, less drop in mean fluid temperature was observed compared to no flow. For the line- and Ltype arrays, clear differences in temperature drop were obtained when considering the flow direction, even for the same magnitude of flow rate. However, in the case of rectangular BHE arrays, the effect of flow direction was almost negligible as expected. The maximum drops in fluid temperature as a function of flow rate were examined in both the least and the most influential flow directions, as shown in Fig. 8. The least influential flow pertains to the directions which have the least affect on the fluid temperature drop (i.e. flow direction angle of 90, 135, and 45 for the line-type, L-type, and rectangular arrays, respectively); the most influential flow pertains to the flow directions which have the minimum reduction in the fluid temperature (i.e. flow direction angle of 0, 45, and 90 for the line-type, L-type, and rectangular arrays, respectively). When the groundwater flow rate was relatively slow or zero, the rectangular BHE array shows a larger drop in fluid temperature than the other types. However, at faster flow rates, it can be observed that the fluid temperature of rectangular type arrays rapidly increased and the difference of fluid temperature due to the array type became smaller, and then became negligible at high flow rates. In terms of sensitivity to the flow rate, the rectangular type array thus seemed to be the most sensitive to the groundwater flow. Regarding to the flow direction, the single linetype array was the most sensitive to the flow direction. When the flow rate was 10 m/yr, the single line-type array gave the largest drop in fluid temperature, as shown in Fig. 7a. The results in Figs. 7 and 8 show that groundwater flow causes less drop in fluid temperature than no flow, which is generally thought to enhance the performance of ground source heat pump systems and are consistent with the results obtained in other researches [18,19,23]. However, at low flow rates (less than 1 m/yr), there are slightly inconsistent results that show almost the same drop in fluid temperature with that of no groundwater flow. To examine this inconsistency, the ground temperature distribution around the BHE arrays was examined. The ground temperature was simulated in the direction of the least influential flow (i.e. the cases shown in Fig. 8a) and the flow rates were selected to be 0.0 (no flow), 0.5 (slow flow), and 10.0 (fast flow) m/yr. When there was no groundwater flow (Fig. 9), it is observed that the rectangular type array resulted in more cold accumulation (more dark area within the contour line of DT ¼ 4.0  C in Fig. 9c than the others) in the ground between boreholes than that for non-rectangular type arrays. That the largest drop in fluid temperature occurs in the rectangular array (see Fig. 8) can be explained by the largest cold accumulation occurring in the rectangular type array, between boreholes. The ground temperature distributions for slow groundwater flow (0.5 m/yr) are shown in Fig. 10. Although there was

236

a

-5.0

Case 1 (Single line type) -5.5 -6.0 -6.5 -7.0 -7.5 -8.0

0

-8.5

0 m/yr (no flow) 0.5 m/yr 1 m/yr 2 m/yr 5 m/yr 10 m/yr

Flow direction Flow direction angle

-9.0 BHEs array

-9.5 0

15

30

45

60

75

90

105

120

135

150

165

180

Maximum drop in mean fluid temperature (° C)

o

Maximum drop in mean fluid temperature ( C)

a

J.C. Choi et al. / Renewable Energy 52 (2013) 230e240

-5.0

least influential flow direction -6.0

-7.0

-8.0

-9.0

-10.0 0.0

o

Flow direction angle ( )

b Case 2 (L -type) -5.5 -6.0 -6.5 -7.0 -7.5 -8.0

0 -8.5

0 m/yr (no flow) 0.5 m/yr 1 m/yr 2 m/yr 5 m/yr 10 m/yr

Flow direction Flow direction angle

-9.0 BHEs array

-9.5 0

15

30

45

60

75

90

105

120

135

150

165

180

o

c

Flow direction angle ( )

Maximum drop in mean fluid temperature (° C)

o

1.0

2.0

3.0

4.0

5.0

6.0

7.0

8.0

9.0 10.0

Groundwater flow rate (m/yr)

b

-5.0

Maximum drop in mean fluid temperature ( C)

o

Case 1 (Single line) θ = 90 o Case 2 (L-type) θ = 135 o Case 3 (Rectangular) θ = 45

-5.0

most influential flow direction -6.0

-7.0

-8.0

o

Case 1 (Single line) θ = 0 o Case 2 (L-type) θ = 45 o Case 3 (Rectangular) θ = 90

-9.0

-10.0 0.0

1.0

2.0

3.0

4.0

5.0

6.0

7.0

8.0

9.0 10.0

Groundwater flow rate (m/yr)

o

Maximum drop in mean fluid temperature ( C)

-5.0

Case 3 (Rectangular type)

Fig. 8. The maximum drop in mean fluid temperature in a) the least and b) the most influential flow direction.

-5.5 -6.0 -6.5

0

-7.0

Flow direction

0 m/yr (no flow) 0.5 m/yr 1 m/yr 2 m/yr 5 m/yr 10 m/yr

Flow direction angle

-7.5 -8.0

BHEs array

-8.5 -9.0 -9.5 0

15

30

45

60

75

90

105

120

135

150

165

180

o

Flow direction angle ( ) Fig. 7. Maximum drop in mean fluid temperature after a 15-yr simulation for a) single line-type, b) L-type, and c) rectangular borehole configurations.

groundwater flow, the ground temperature distributions were similar to that without groundwater flowing, which are shown in Fig. 9. This may explain why very low flow rates resulted in a negligible change of fluid temperature (see Fig. 8). Similar results were obtained by Signorelli et al. [21], who found negligible effects of relatively low groundwater flow rates on effective thermal conductivities using numerical evaluation of TRT. To quantify the

effect of heat advection vs. conduction, Peclet number, Pe ¼ urwcwLc/lg, corresponding to its flow rate was examined. The Peclet number can be a good indicator to quantify the effect of ground advection. However, it can vary with the characteristic length Lc. Chiasson et al. [16], Sutton et al. [17], Diao et al. [18], and Molina-Giraldo et al. [36] assumed the characteristic length as the spacing between the boreholes, the borehole radius, the radial distance, and the height of the borehole, respectively. For consistency, the characteristic length of Pe in this study was assumed to be a unit number (i.e. Lc ¼ 1). Calculated Pe corresponding to the very low flow velocity (less than 1 m/yr) was around 0.05. When the Pe in the aforementioned research [16e18,36], which define the criteria of negligible advection effect, were re-calculated as Pe with unit characteristic length, the Pe less than 0.05 which was obtained by this study was within the range of previously reported values (0.02e0.22). This result indicates that, regardless of the array type, low rates of groundwater flow (less than around 1 m/yr or Pe is less than 0.05) appear not to affect heat transfer around heat exchangers, and thus conduction would be the prevailing heat transfer mechanism. The ground temperature distributions generated using a flow rate of 10 m/yr are shown in Fig. 11. The cold fronts for the three cases, where DT ¼ 0.5  C for the ground temperature, were

J.C. Choi et al. / Renewable Energy 52 (2013) 230e240

Fig. 9. The distribution of maximum ground temperature change around heat exchangers with no groundwater flow for a) single line-type, b) L-type, and c) rectangular borehole arrays.

237

Fig. 10. The distribution of maximum ground temperature change around heat exchangers with slow groundwater flow (0.5 m/yr) for a) single line-type, b) L-type, and c) rectangular borehole arrays.

238

J.C. Choi et al. / Renewable Energy 52 (2013) 230e240

Fig. 11. The distribution of maximum ground temperature change around heat exchangers with fast groundwater flow (10.0 m/yr) for a) single line-type, b) L-type, and c) rectangular borehole arrays.

J.C. Choi et al. / Renewable Energy 52 (2013) 230e240

350

Least influential direction Most influential direction

Annual heating capacity (MWh)

300

13%

6%

250 200 150 100 50 0 Single-line

L-type

Rectangular-type

Fig. 12. Comparison of annual heating capacity associated with flow direction when the flow rate is 10 m/yr.

similar. However, when local cold accumulations between boreholes were compared, the single line array resulted in a larger change in the ground temperature near the downstream boreholes (see the larger gray contour of DT ¼ 2.0  C in Fig. 11a) compared to the other arrays. This larger amount of local cold accumulation for the single line array seemed to be related to the number of boreholes parallel to groundwater flow. The opposite results can be observed in Fig. 8b, where the single line array gave the smallest drop in fluid temperature when each borehole was perpendicular to the flow. Diao et al. [18] reported a similar result in a simulation with a 2  3 borehole array, where a single line of boreholes perpendicular to the flow direction was found to be the most effective to enhance the performance. In case of the ground thermal storage system (e.g., boreholes coupled with solar system to inject the heat into the ground), more heat/cold accumulation can be beneficial. However, in general ground heat exchanger systems, minimizing the number of boreholes along the flow direction should be ensured to enhance the performance of a ground source heat pump system. The annual heating capacity associated with flow direction was compared as shown in Fig. 12. A seasonal performance factor (SPF) for heating was selected to be a value of 3.0. Monthly distribution ratios of total and ground energy loads were assumed to be the same as that in Fig. 3. The constraint for the minimum fluid temperature was set to be 0  C. The annual heat capacity considering the most influential flow direction for line-type, L-type, and rectangular borehole arrays were 295, 291, and 275 MWh, respectively. Those considering the least influential flow direction were 262, 275, and 275 for line-type, L-type, and rectangular borehole arrays, respectively. The annual heating capacity was up to 13% greater when comparing the most influential flow direction with the least influential flow. This implies that consideration of the direction of groundwater flow can save around 10% of the installation and operating costs. 5. Conclusion The influence of both direction and rate of groundwater flow on different types of BHE arrays were investigated numerically. The coupled heat conductioneadvection model based on the FE method was adopted for the heat simulation outside boreholes, and the heat transfer inside the boreholes was treated as steady-state thermal resistance. Long-term performances considering 15-yr heating operations were simulated on three different types of BHE arrays: single line, L-type, and rectangular. The model was then

239

compared with the moving infinite line source and the finite line source model and gave good agreement. Non-square rectangular arrays were noticeably influenced by the direction of groundwater flow. The single line array was mostly influenced by the direction of the groundwater flow. In terms of sensitivity to the flow rate (or the magnitude of flow velocity), the rectangular type array seemed to be the most sensitive, compared to the other arrays. The comparison of the annual heat capacity showed differences up to 13% when considering the flow direction. This result indicates that the consideration of the direction of groundwater flow to estimate the size of BHE arrays may be important, as well as considering flow rates, particularly for the non-square rectangular array. When the Peclet number (with the unit characteristic length) is less than 0.05, the temperature distribution in the vicinity of borehole arrays shows that groundwater flows seem not to affect the temperature distribution around boreholes, regardless of array type. However when the Peclet number is higher than 0.05 (i.e. the flow rate is higher than 1 m/yr in this study), the groundwater flow could induce local cold accumulations around some boreholes, particularly in the downstream direction, and also seemed to enhance the drop in fluid temperature. When the positive effects of groundwater flow on a non-axisymmetric array is intended, both low groundwater flow rate and numbers of heat exchangers along the flow direction should be avoided. Although the importance of groundwater flow direction on the various types of BHE arrays was roughly examined numerically in this study, the heat transfer in the longitudinal direction of the borehole axis could affect the performance of the system for longterm operation and the transient heat transfer inside boreholes may affect the short-term performance of the BHE arrays. Future work should therefore include effects of the finite length of the borehole and the transient heat transfer inside boreholes in the numerical model. In addition, further integrated research based on field experimental approaches with more realistic energy load (e.g., peak heating demands, non-constant heat flux along the vertical dimension and each borehole, etc.) and geological condition (e.g., in-homogeneity of the ground) may give better understanding of the effects of groundwater flow on the BHE arrays. Acknowledgment This study was supported by the Research Council of Norway through the NGI (Norwegian Geotechnical Institute)’s postdoctoral research fellowship and the 2011 Construction Technology Innovation Project (11E04) by the KICTEP under the Ministry of Land, Transport, and Maritime Affair of Korea. References [1] Lund JW. Direct heat utilization of geothermal resources. Renew Energ 1997; 10:403e8. [2] Sanner B, Karytsas C, Mendrinos D, Rybach L. Current status of ground source heat pumps and underground thermal energy storage in Europe. Geothermics 2003;32:579e88. [3] Yang W, Zhou J, Xu W, Zhang G. Current status of ground-source heat pumps in China. Energ Policy 2010;38:323e32. [4] Midttømme K, Berre I, Hauge A, Musæus TE, Kristjànsson BR. Geothermal energy e country update for Norway. In: Proceedings of the 2010 World Geothermal Congress; 2010. p. 25e9. [5] Ingersoll LR, Zobel OJ, Ingersoll AC. Heat conduction with engineering, geological and other applications. New York: Mcgraw-Hill; 1954. [6] Carslaw HS, Jaeger JC. Conduction of heat in solids. Claremore Press; 1959. [7] Eskilson P. Thermal analysis of heat extraction boreholes. Department of Mathematical Physics, University of Lund; 1987. [8] Zeng H, Diao N, Fang Z. Heat transfer analysis of boreholes in vertical ground heat exchangers. Int J Heat Mass Transf 2003;46:4467e81. [9] Lamarche L, Beauchamp B. A new contribution to the finite line-source model for geothermal boreholes. Energ Buildings 2007;39:188e98.

240

J.C. Choi et al. / Renewable Energy 52 (2013) 230e240

[10] Hellström G. Ground heat storage thermal analysis of duct storage systems. Part I. Theory. Department of Mathematical Physics, University of Lund; 1991. [11] Yavuzturk C, Spitler JD. A short time step response factor model for vertical ground loop heat exchangers. ASHRAE Trans 1999;105:475e85. [12] Muraya NK, O’Neal DL, Heffington WM. Thermal interference of adjacent legs in a vertical U-tube heat exchanger for a ground-coupled heat pump. ASHRAE Trans 1996;102:12e21. [13] Li Z, Zheng M. Development of a numerical model for the simulation of vertical U-tube ground heat exchangers. Appl Therm Eng 2009;29:920e4. [14] Choi JC, Lee SR, Lee DS. Numerical simulation of vertical ground heat exchangers: intermittent operation in unsaturated soil conditions. Comput Geotech 2011;38:949e58. [15] Witte H. Geothermal response tests with heat extraction and heat injection: examples of application in research and design of geothermal ground heat exchangers. In: Workshop geothermal response tests. Lausanne: Ecole Polytechnique Féd’erale Lausanne, Schweizerische Vereinigung Für Geothermie SVG; 2001. p. 48e63. [16] Chiasson AD, Rees SJ, Spitler JD. A preliminary assessment of the effects of ground-water flow on closed-loop ground-source heat pumps systems. ASHRAE Trans 2000;106:380e93. [17] Sutton MG, Nutter DW, Couvillion RJ. A ground resistance for vertical bore heat exchangers with groundwater flow. J Energ Resour e ASME 2003;125: 183e9. [18] Diao N, Li Q, Fang Z. Heat transfer in ground heat exchangers with groundwater advection. Int J Therm Sci 2004;43:1203e11. [19] Fan R, Jiang Y, Yao Y, Shiming D, Ma Z. A study on the performance of a geothermal heat exchanger under coupled heat conduction and groundwater advection. Energy 2007;32:2199e209. [20] Gehlin SEA, Hellström G. Influence on thermal response test by groundwater flow in vertical fractures in hard rock. Renew Energ 2003;28:2221e38. [21] Signorelli S, Bassetti S, Pahud D, Kohl T. Numerical evaluation of thermal response tests. Geothermics 2007;36:141e66. [22] Raymond J, Therrien R, Gosselin L, Lefebvre R. Numerical analysis of thermal response tests with a groundwater flow and heat transfer model. Renew Energ 2011;36:315e24. [23] Fujii H, Itoi R, Fujii J, Uchida Y. Optimizing the design of large-scale groundcoupled heat pump systems using groundwater and heat transport modeling. Geothermics 2005;34:347e64. [24] Yang W, Shi M, Liu G, Chen Z. A two-region simulation model of vertical Utube ground heat exchanger and its experimental verification. Appl Energ 2009;86:2005e12.

[25] Kim EJ, Roux JJ, Rusaouen G, Kuznik F. Numerical modelling of geothermal vertical heat exchangers for the short time analysis using the state model size reduction technique. Appl Therm Eng 2010;30:706e14. [26] Al-Khoury R, Kölbel T, Schramedei R. Efficient numerical modeling of borehole heat exchangers. Comput Geosci 2010;36:1301e15. [27] Woodside W, Messmer JH. Thermal conductivity of porous media. I. Unconsolidated sands. J Appl Phys 1961;32:1688e99. [28] Beck AE. An improved method of computing the thermal conductivity of fluidfilled sedimentary rocks. Geophysics 1976;41:133e44. [29] Marcotte D, Pasquier P. On the estimation of thermal resistance in borehole thermal conductivity test. Renew Energ 2008;33:2407e15. [30] Bose J, Parker JD, McQuiston FC. Design/data manual for closed-loop groundcoupled heat pump systems. Atlanta: American Society of Heating, Refrigerating, and Air-Conditioning Engineers; 1985. [31] Remund CP. Borehole thermal resistance: laboratory and field studies. ASHRAE Trans 1999;105:439e45. [32] Lamarche L, Kajl S, Beauchamp B. A review of methods to evaluate borehole thermal resistances in geothermal heat-pump systems. Geothermics 2010;39: 187e200. [33] Incropera F, DeWitt D. Fundamentals of heat and mass transfer. New York: Wiley; 1996. [34] Zienkiewicz OC, Taylor RL, Nithiarasu P. The finite element method for fluid dynamics. ButterwortheHeinemann; 2005. [35] EED 3.0 manual e Earth energy designer. Blocon; 2010. [36] Molina-Giraldo N, Blum P, Zhu K, Bayer P, Fang Z. A moving finite line source model to simulate borehole heat exchangers with groundwater advection. Int J Therm Sci 2011;50:2506e13. [37] Rybach L, Sanner B. Ground source heat pump systems, the European experience. Geo-Heat Cent Q Bull 2000:16e26. [38] Lee CK, Lam HN. Computer simulation of borehole ground heat exchangers for geothermal heat pump systems. Renew Energ 2008;33:1286e96. [39] Philippe M, Bernier M, Marchio D. Validity ranges of three analytical solutions to heat transfer in the vicinity of single boreholes. Geothermics 2009;38:407e13. [40] Marcotte D, Pasquier P, Sheriff F, Bernier M. The importance of axial effects for borehole design of geothermal heat-pump systems. Renew Energ 2010;35: 763e70. [41] Domenico PA, Schwartz FW. Physical and chemical hydrogeology. New York: Wiley; 1990. [42] Cook P. A guide to regional groundwater flow in fractured rock aquifers. Australia: CSIRO Land and Water; 2003.