Performance of heat pumps with direct expansion in vertical ground heat exchangers in heating mode

Energy Conversion and Management 95 (2015) 120–130

Contents lists available at ScienceDirect

Energy Conversion and Management journal homepage: www.elsevier.com/locate/enconman

Performance of heat pumps with direct expansion in vertical ground heat exchangers in heating mode Michele De Carli ⇑, Stefano Fiorenzato, Angelo Zarrella Department of Industrial Engineering, Applied Physics Section, University of Padova, Via Venezia 1, 35131, Italy

a r t i c l e

i n f o

Article history: Received 7 October 2014 Accepted 29 January 2015

Keywords: Ground source heat pumps Direct expansion Convective heat transfer coefficient Borehole heat exchanger Coefficient of performance

a b s t r a c t Ground source heat pump systems represent an interesting example of renewable energy technology for heating and cooling of buildings. The connection with the ground is usually done by means of a closed loop where a heat-carrier fluid (pure water or a solution of antifreeze and water) flows and, in heating mode, moves heat from ground to refrigerant fluid of heat pump. A new solution is the direct expansion heat pump. In this case, the heat-carrier fluid inside the ground loop is the same refrigerant fluid of heat pump. This paper focuses on the energy performance of direct expansion ground source heat pump with borehole heat exchangers in heating mode, looking at residential building installations. For this purpose, the evaporating process of the refrigerant fluid inside vertical tubes is investigated in order to analyze the influence of the convective heat transfer coefficient on the global heat transfer with the surrounding ground. Then, an analytical model reported in literature for the design of common borehole heat exchangers has been modified for direct expansion systems. Finally, the direct expansion and common ground source heat pumps have been compared in terms of both total borehole length and thermal performance. Results indicate that the direct expansion system has higher energy performance and requires lower total borehole length compared to the common system. However, when the two systems are compared with the same mean fluid evaporating temperature, the overall length of the ground heat exchanger of the direct expansion heat pump is greater than that of the common system. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction High energy efficiency and decrease of emissions are essential objectives for the design of space conditioning systems. Ground source heat pumps are one of the most suitable solutions for this purpose, due to their high energy performance, named COP (i.e. Coefficient of Performance) [1–3]. In this case, soil is used as a ‘‘hot sink’’ or ‘‘cold sink’’, or rather, a thermal source where heat is absorbed or injected. The advantage is that the ground temperature is more constant and uniform than ambient air temperature [4]. The connection between the heat pump and soil is possible for example with closed loops, named Secondary Loop Ground Source Heat Pumps (SL-GSHPs). The design and sizing of these systems is always debated and several methods are widely discussed in literature. Recently attention has been addressed to systems with a direct connection between the heat pump and the soil, namely Direct Expansion Ground Source Heat Pumps (DX-GSHPs). For these ⇑ Corresponding author. Tel.: +39 049 827 6870; fax: +39 049 827 6896. E-mail address: [email protected] (M. De Carli). http://dx.doi.org/10.1016/j.enconman.2015.01.080 0196-8904/Ó 2015 Elsevier Ltd. All rights reserved.

systems, the design of ground heat exchangers (GHEs) is delicate and difficult to analyze, due to the physical phase transition of the refrigerant fluid inside buried pipes. There are few experimental studies about DX-GSHP. Wang et al. [5] conducted experimental tests in heating mode with a direct expansion ground source heat pump which used R134a as refrigerant; the COP values of the heat pump and the entire system were found to be on average 3.55 and 3.28, respectively, when the evaporating temperature was 3.14 °C and condensing temperature was 53.4 °C. Wang et al. [6] analyzed a direct expansion ground source heat pump with R22 as refrigerant in heating mode and installed in China; they found an average COP of the heat pump equal to 3.12 and a heat-absorption rate per unit of GHE length of about 54.4 W/m. Fannou et al. [7] presented an experimental analysis of a direct expansion ground source heat pump with R22 as refrigerant installed in Montreal (Canada); their test measurements were conducted over a one-month in early spring and showed COP values between 2.7 and 3.44. Beauchamp et al. [8] proposed a model to analyze the transient behavior of a vertical ground heat exchanger by an appropriate thermal resistances and capacities system. The U-tube is discretized in control volume and, for each element, mass and energy

M. De Carli et al. / Energy Conversion and Management 95 (2015) 120–130

121

Nomenclature cp dbore did dod dmax f Fo G GF GHE k L _ m n N P Pr Q_ r R Rb Rga Rgd Rgm Re s T Tsoil x

specific heat (J kg1 K1) borehole diameter (m) pipe inner diameter (m) pipe outer diameter (m) diameter from borehole axis beyond which the undisturbed ground temperature is assumed (m) friction factor Fourier number mass flow rate of fluid per cross-sectional area of flow (kg m2 s1) G-factor, used for determining soil thermal resistance ground heat exchanger thermal conductivity (W m1 K1) length (m) fluid mass flow rate (kg s1) number of boreholes number of steps used for evaporating process power (W) Prandtl number thermal power (W) specific latent heat (J kg1) thermal resistance per unit length (m K W1) filling thermal resistance per unit length (m K W1) thermal resistance of the ground related to annual impulse per unit length (m K W1) thermal resistance of the ground related to daily impulse per unit length (m K W1) thermal resistance of the ground related to monthly impulse per unit length (m K W1) Reynolds number thickness (m) temperature (K) undisturbed soil temperature (K) vapor quality

balance equations were written; as a matter of fact, the model solves the transient problem using a purely implicit finite difference method. The results were temperature profiles, for both flowing refrigerant fluid and pipe wall, considering also the transient thermal interference between the pipes of the U-loop. The thermodynamic process, and hence, the convective heat transfer coefficient, as well as temperature and pressure losses were not evaluated because of the complexity of the two-phase flow related to this approach. Austin and Sumathy [9] presented a parametric study on the performance of a direct-expansion ground source heat pump using carbon dioxide and horizontal buried pipe. In their approach each component of the machine was modelled and the GHE was the evaporator of the heat pump; thermal resistances were used to simulate the heat exchange with the ground and, at the same time, the convective heat transfer coefficient and the pressure losses were estimated. Using that approach, it was possible to link the heat exchange and the thermodynamic process of the refrigerant but dynamic behavior and thermal interference between pipes were not taken into account. Eslami-Nejad et al. [10] developed a numerical model to evaluate the thermal performance of a carbon dioxide filled vertical ground heat exchanger consisting of a long copper U-tube; their study did not include the heat pump. According to their results, large amounts of energy were absorbed from the ground due to the high two-phase heat transfer coefficient of CO2. The drill’s analysis has two main aspects that may be enforced in the model: the first one is an accurate characterization of the

Greek symbols a convective heat transfer coefficient (W m2 K1) e vapor volume fraction uLO two-phase multiplier (–) g efficiency (–) q density (kg m3) r surface tension (Pa m2) l dynamic viscosity (kg m1 s1) Subscripts c critical point comp compression eb boiling point el electrical ev evaporating process fr friction g ground, system of ground heat exchangers gr gravity hD heating design condition i index corresponding to a generic state of the process, evaporation or superheating L liquid, phase fluid considered flows only with the own mass flow rate LO liquid phase flows with the total mass flow rate mech mechanical mix liquid–vapor mixture ov superheating process p pipe wall R refrigerant SAT saturation point sys system V vapor, phase fluid considered flows only with the own mass flow rate

heat exchange with thermal behaviour and interference phenomena within the drills; the second one is the thermodynamic process of the refrigerant fluid, with the calculation of the fluid convective heat transfer coefficient, the pressure and the temperature trends for steady state conditions. The design of GHEs is usually similar for both the traditional systems (SL-GSHP) and the DX-GSHP: it is possible to follow and adapt the same approach to model soil and borehole, while the refrigerant has to be considered in a suitable way. The basic models use analytical solution of line source [11] and cylindrical source [12]. Using these approaches, Kavanaugh and Rafferty [13] proposed a design method of vertical GHEs (i.e. the well-known ASHRAE method), which considers the thermal interference between adjacent boreholes by means of a ‘‘ground penalty temperature’’. Eskilson [14] proposed dimensionless parameters, named g-functions, to describe the performance of boreholes in homogenous ground. Hellström and Sanner [15] presented another approach which combines analytical and numerical methods to calculate the borehole thermal resistance and uses the solution of cylindrical source given by Carslaw and Jaeger [12] to estimate thermal resistance between the borehole wall and the undisturbed ground. These models are valid when the long-term behavior is analyzed, i.e. in order to investigate the effect of the thermal drift of the ground or the interaction between boreholes [16]. When the short-time behavior (between 2 and 6 h) is considered, the previous approaches may give wrong results since the borehole thermal capacity is not modelled [17,18].

122

M. De Carli et al. / Energy Conversion and Management 95 (2015) 120–130

This present study shows an analysis of the ground heat exchanger in the long-term for the heating season. The model considers in detail all aspects of the thermodynamic process of the two-phase problem, combined with the calculation of the heat exchange with the ground. As will be shown afterward, first a steady state model has been carried out in order to check the coupling of the two-phase evaporating problem inside the pipes with the model of the ground heat transfer. Then, a more detailed model has been developed using an appropriate thermal resistance system based on the method to design ground vertical boreholes proposed by Kavanaugh and Rafferty [13]. As a final step, the heat pump cycle has been integrated in the overall model, thus allowing to simulate the entire DX-GSHP system, which has been compared to a SL-GSHP.

Pipe

Entering fluid (i)

Grout Rb

Soil

Undisturbed conditions

2. Theory For describing the complete model it is useful to subdivide the problem in three main tasks, as reported in detail hereafter:  The model for the design of the borehole (from the internal surface of the pipes to the undisturbed ground through the filling material of the borehole).  The model of the refrigerant fluid to determinate the convective heat transfer coefficient on the internal wall of the pipes.  The model of the refrigerant fluid to determine the pressure drops along the pipes. The present study looks at the long-term performance of DXGSHP systems during the heating season. In particular, the main goal is the design of the DX-GHE coupled with the heat pump, evaluating at the same time the thermodynamic process inside the pipe under the following initial hypotheses:  single U-tube borehole heat exchanger;  heat transfer in the soil occurs only by conduction and ground water flow is neglected;  only one equivalent ground layer is considered;  filling material thermal conductivity is uniform and constant;  soil vertical temperature profile is uniform and equal to the undisturbed temperature;  thermal capacity of the refrigerant fluid is neglected;  as mentioned in the introduction, soil thermal capacity has been in a first step neglected in order to check the coupling of the models of the refrigerant fluid and GHE; in the second phase of the work it has been included, according to [13];  in the first stage of the analysis, the thermal interference between the supply (downward) and return (upstream) pipes has been neglected, as well as possible interferences within the borehole field; in the second part of the work they have been included.

Outgoing fluid (i+1) Fig. 1. Scheme of the thermal model of the ith segment describing the heat exchange between the refrigerant fluid within the pipe and the ground in the simplified model.

ðT soil  T R;i Þ  ðT soil  T R;iþ1 Þ Q_ ði;iþ1Þsys ¼ K  Ai;iþ1  Li;iþ1  ðT T R;i Þ ln ðT soilT R;iþ1 Þ soil

¼ Q_ ði;iþ1ÞR

ð1Þ

where

KAi;iþ1 ¼

ln ddod ln ddmax 1 1 id bore þ þ þ Rb 2p  did  ai R 4p  kp p  ksoil 2

!1 ð2Þ

Rb is the thermal resistance of the grouting material, calculated making use of the relationship proposed by Paul and Remund [19,20], while Q_ sys represents the heat flux rate exchanged by the overall system refrigerant-borehole-soil. Along the pipe length two regions are present, i.e. the two-phase region and the superheated vapor region. Based on this consideration, Q_ iðRÞ may be expressed by means of Eq. (3) for the two phase region and Eq. (4) for the superheated vapor region:

_  rðT R Þ  Dxði;iþ1ÞR Q_ ði;iþ1ÞR ¼ m

ð3Þ

_  cp  DT ði;iþ1ÞR Q_ ði;iþ1ÞR ¼ m

ð4Þ

where r(TR) is the specific latent heat at the temperature TR, which can be calculated with the following equation [21]:



Tc  TR T c  T eb

0:38

2.1. Simplified borehole model

rðT R Þ ¼ r eb 

In the present section the first stage of the work is presented, i.e. a simplified model for describing the heat exchange between the U-tube, the grout material and the ground has been considered (neither thermal short-circuit between the pipes inside the bore nor thermal capacity of the ground). Under the above mentioned hypotheses, considering the U-tube as only one pipe subdivided into segments along the axis (Fig. 1), between the generic inlet and outlet fluid sections, it is possible to write for the ith segment the following equation:

This model has to be linked to the evaluation of the convective heat transfer coefficient and the pressure losses related to the effective refrigerant thermodynamic process.

ð5Þ

2.2. Refrigerant fluid: convective heat transfer coefficient Inside the evaporator, two regions are present and the evaluation of the convective heat transfer coefficient depends on the type of refrigerant flow realized inside the heat exchanger, i.e. this

123

M. De Carli et al. / Energy Conversion and Management 95 (2015) 120–130

coefficient depends on the motion type inside the pipe (i.e. annular, bag or spray flows) [22]. Considering, for example, a vertical pipe subject to a heat flux in the evaporating process and fed by a constant mass flow rate, different kinds of two-phase outflow layouts may happen [22]. If the refrigerant enters as subcooled liquid, hence heat exchange happens following the forced convection within the liquid phase, the heat transfer coefficient is constant, assuming unchanged thermal properties with the temperature; in this case temperature profiles are linear along the length of the pipe and the difference between pipe wall temperature and refrigerant temperature is constant. At any distance from the inlet section, pipe wall temperature allows the bubble’s formation and the phenomenon is known as nucleate boiling within subcooled liquid. In this region pipe wall temperature is constant and the difference between the wall and the refrigerant temperatures decreases if the distance increases. Until boiling is related to saturated liquid, bubble and sac outflows are present, leading to a nearly constant convective heat transfer coefficient. As the evaporating process evolves, refrigerant flow becomes annular and a new type of boiling phenomenon starts: heat transfer happens due to the forced convection within the liquid film up to the vapor-liquid interface (‘‘forced two-phase boiling’’). As the liquid film thickness decreases, nucleation sites disappear and the boiling is suppressed gradually until the complete evaporation. In this work, the complex characterization of the two-phase outflow has been defined via Chen’s correlation [23,24], since it covers the nucleate boiling liquid saturated region and the forced twophase boiling region. According to this theory, the convective heat transfer coefficient is calculated as a sum of the nucleate boiling contribute and the forced two-phase convection contribute, assessed with appropriate coefficients that represent the trend of the process: 

aR ¼ 0:023  " þ



kL  ðReL Þ0:8  ðPrL Þ0:4  F did 0:79

0:00122  ðkL Þ

 ðcpL Þ0:45  ðqL Þ0:49  ðT p  T SAT Þ0:24  ðDpSAT Þ0:75

ðrÞ

0:5

 ðlL Þ

0:29

 ðr SAT Þ

0:24

 ðqV Þ

0:24

if

1=X tt  0:1

F ¼ 2:35  ð1=X tt þ 0:213Þ0:736

if

1=X tt > 0:1

# S

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi u dp " #0:5 u u dz fr;L f L ð1  xÞ2 qV u   X tt ¼ t ¼   fV x2 qL  dp dz

ð6Þ

ð7Þ

ð8Þ

fr;V

The coefficient S is a suppression factor to account for the reduction in the coefficient due to the suppression of nucleation sites; it is:

S¼

1 1 þ 2:53  106  ðReL  F 1:25 Þ

1:17

ð9Þ

In the superheated vapor region the convective heat transfer coefficient is calculated by means of the relationship proposed by Gnielinski and Petukhov [24]:

aR ¼

f 8

 ðReV  1000Þ  PrV kv   0:5 did 2=3 f 1 þ 12:7  8  ½ðPrÞ  1

ð10Þ

with the friction factor f given by the following equation:

f ¼ ½0:79  lnðReV Þ  1:64

2

Once the convective heat transfer coefficient is estimated, pressure losses have to be considered since they involved a decrease of the saturation temperature, consequently variations of the thermal properties of the refrigerant occur. The pressure losses depend on the adopted liquid–vapor layout. Two models can be considered: the homogenous outflow model or the separated-phase outflow model. The last model is more generic and accurate, consequently allows to simplify the problem [21]. For this purpose, it is necessary to value the ‘‘volume fraction occupied from vapor’’ (e), calculating the gas-phase velocity and liquid-phase velocity ratio and considering adequately the mass flow rate [22]. Pressure losses may be distinguished into three components: friction, gravity and amount of the motion variation component. They are calculated considering viscosity and density referred to the refrigerant temperature ‘‘TR’’ at the ‘‘i’’ quality. As for friction losses, in the two-phase region they are due only to liquid phase because of the annular motion inside the pipe.

ðDpi Þfr ¼ 2  f LO 

Li;iþ1 G   ðULO Þ2 did qL

ð12Þ

where fLO is the liquid-only friction factor which can be calculated by Eq. (13) [25] and ULO is the liquid-only two-phase multiplier calculated using the correlation proposed by McAdams et al. [26] as indicated in Eq. (14):

f LO ¼

0:0791 ðReL Þ0:25

ReL  2300

ð13Þ

       0:1 lL qV ðULO Þ2 ¼ 1 þ  1  xi  1 þ  1  xi

lV

qL

ð14Þ

Gravity losses are calculated based on the integral in the considered process length Li:

The coefficient F represents the two-phase flow based on the Martinelli’s parameter Xtt:

F¼1

2.3. Refrigerant fluid: pressure losses

ð11Þ

ðDpi Þgr ¼

Z

iþ1

ðg  qmix Þdz

i

ð15Þ

with

qmix ¼ e  qV þ ð1  eÞ  qL

ð16Þ

The pressure losses due to the motion variation are based on the integral in the considered process length Li:

ðDpi Þa ¼ G2 

Z i

iþ1



 x2 ð1  x2 Þ dz þ e  qV ð1  eÞ  qL

ð17Þ

For the superheated vapor region the over-written equations can be applied by setting as 1 both the vapor quality and the volume fraction occupied by the vapor. 2.4. Improved model for GHE The design of a GSHP requires the design peak load, thermal energy needs of the building and energy efficiency (i.e. COP) of the heat pump. In this way, based on the physical parameters of the soil and once decided the geometry and the characteristics of the borehole, the overall length of the ground heat exchangers can be determined. Usually the length depends on both heating and cooling conditions, but, as already mentioned, in the present work the model considers only the evaporating process of the refrigerant. The improved method for analyzing the GHE thermal behavior (Fig. 2) is based on the analytical method [13], which has been coupled to the refrigerant evaporation model, as described hereafter. As for the mass flow rate, once fixed the number of buried drills n and the inlet vapor quality, the following equation can be written:

124

M. De Carli et al. / Energy Conversion and Management 95 (2015) 120–130

Downward stream

  Dp _  cpV  ðT cond  T R;out Þ þ n  m _   Picomp ¼ n  m

Interference between pipes

qV

1

gel  gmech

ð22Þ

The whole inputs and outputs of the heat pump are known, therefore the COP of the system can be determined once calculated the sizing of the GHEs. The COP can be then substituted at the beginning of the program in a iterative process. 3. Methods

Upward stream

Interference between boreholes

Fig. 2. Scheme of the thermal model of two parallel generic segments describing the heat exchange between the refrigerant fluid within the pipe and the ground in the improved model.

_ ¼ m

Q_ g;hD n  rðT R Þin  Dx

ð18Þ

where Q_ g;hD is the total heat load exchanged with the ground in design conditions. The analytical method for the ground heat exchange is based on the annual heat demand of the building Q_ a , the partial load factor (PLFhD) and the thermal short-circuit factor (Fsc) [13]. The approach is based on the soil response to three appropriate temporal impulses (named G-factors), which are evaluated by means of the Fourier number and the following equation:

GF ¼ 0:0758  lnðFoÞ þ 0:1009

ð19Þ

Consequently, the thermal resistance of the ground related to annual (Rga), monthly (Rgm) and daily (Rgd) impulses are determined, while the overall sum of the thermal resistances (except Rga) is expressed as follows:

X

  ln ddod 1 Rb id þ þ þ 2  PLF hD  Rgm þ 2  F sc  Rgd Ri;iþ1 ¼ 2  p  did  aiR 4  p  kp 2 ð20Þ

The expression of the required length of the GHEs has to consider the above thermal resistances, the effect of the annual heat flux ðQ_ a Þ and the evaporating process:

Li;iþ1 ¼

P _  rðT i R Þ  Dxði;iþ1ÞR  Rði;iþ1Þ m Q_ a  Rga þ T i soil  T i R  T pen N  ðT i soil  T i R  T pen Þ

For the simulation of the drill behavior a Fortran code has been developed. Fig. 3 shows the flow chart that resumes the main steps of the program. The convective heat transfer coefficient as well as the pressure and temperature drops are evaluated; the model described in the Section 2.1 is enforced to obtain the heat transfer fluxes and the required length of the pipe (i.e. the size of the evaporator). For this purpose the evaporating process is divided into N steps in order to recognize the evaporating and superheating zones. Consequently, the core of the code is divided into two main sections that simulate the behavior of the two regions existing inside the evaporator. The program requires the following input data: geometrical characteristics and thermal conductivity of the soil and of the borehole, mass flow rate, evaporation initial quality, U-tube total length and number of steps N; the vapor quality is considered an independent variable. Depending on temperature and pressure, refrigerant properties are calculated and updated step by step. In a further step, pressure losses are evaluated as well as the length Li,i+1. At this point it is possible to estimate the new pressure and temperature conditions for the stage (i + 1). In case of complete evaporation (x = 1), superheated vapor is considered. Using this approach, it is possible to estimate the evaporation fluid behavior in the pipe. As for the model of the heat transfer from the pipe to the ground through the grouting, two models have been developed. First, a simplified steady state model of the heat transfer in the ground has been considered in order to check the contemporary calculation of the evaporating process and the heat transfer with the undisturbed ground. Then, a more detailed model for the ground has been developed, which allows to consider the analytic method for sizing the GHE length. Finally, a comprehensive model which considers the GHE and the overall heat pump cycle has been carried out for evaluating the performance of the entire system, as shown afterward. Before using the developed model, the influence of the number of steps N was investigated. Fig. 4 shows the profile of the convective heat transfer coefficient for R410a and different values of the parameter N, maintaining unchanged other input data. As can be seen, the effect of N is negligible; however high values of this parameter allows a better pattern and separation of the two regions through the ground heat exchanger (i.e. two-phase and superheated vapor regions). The same outcome was also found for saturation temperature and pressure. 4. Results

ð21Þ

where Tpen is the ground penalty temperature [13,27] which allows the evaluation of thermal interference between adjacent boreholes. Finally, the evaluation of the thermodynamic process of the refrigerant has been introduced for the simultaneous solution. As the outlet conditions and the condensing temperature are output of the calculation, it is possible to estimate the compressor power. In the present model the heat exchange due to the horizontal piping between GHEs and the heat pump is neglected, while the pressure drop is considered. As for the compression process, ideal gas behavior of the refrigerant is supposed. With these assumptions the power of the compressor is expressed as follows:

4.1. Simulation of the evaporating process The first analysis concerns the profiles trend of the refrigerant convective heat transfer coefficient; according to Chen’s theory [23,24], it depends on the specific mass flow rate, saturation temperature and quality. The first test refers to R410A, inlet vapor quality x = 0.1 and a discretization of the pipe into N = 500 elements; results reported in Fig. 5 show the profiles of the convective heat transfer coefficient, related to the vapor quality, varying inlet temperature of the refrigerant for a mass flow rate of the fluid equal to 0.012 kg s1. The values of the heat transfer coefficient increase

125

M. De Carli et al. / Energy Conversion and Management 95 (2015) 120–130

INPUT data: characteristics borehole-soil, ṁ, x in , TR in ,pR in ,N, Ltot

Initialization quality vector

i=i+1

CHECK: cumulated length and quality if : L ev< Ltot & xi < 1

if : L ev < Ltot & xi = 1

CHECK: cumulated length and refrigerant temperature

procedure : EVAPORATING FLUID

if : L ov < Ltot & TiR < Tmax

Procedure : SUPERHEATING FLUID • • •

•

Convective heat transfer coefficient valuation αiR

•

Pressure losses evaluation (Pa/m) following separated-phase model Length of the process computation

•

Li,i+1 , cumulated length valuation L ev , total heat transfer, next TSAT and pSAT

Convective heat transfer valuation Pressure losses evaluation (Pa/m) Length of the process computation Li,i+1,

for stage (i+1)

cumulated length L ov and total heat transfer .

L ev= Ltot & xi = 1

END PROCESS output profiles : αiR ,

L ov= Ltot

TSAT , pSAT , total heat transfer

END PROCESS output profiles: αiR , T, p, total heat flux

Fig. 3. Scheme of the calculation steps for the refrigerant fluid in the GHE.

2500

convecve heat transfer coefficient N = 100

N = 250

N = 500

inlet temperature -5°C inlet temperature 5°C

N = 700

2500 2000

1500

α [W/(m2 K)]

α (W/(m2K))

2000

1000 500 0

inlet temperature 0°C inlet temperature 10°C

1500 1000 500

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

quality

0 0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

quality

Fig. 4. Effect of the number of steps N on the convective heat transfer coefficient for R410A for a mass flow rate of 0.012 kg/s.

Fig. 5. Profiles of the convective heat transfer coefficient for R410A at different inlet temperatures for the whole evaporating process for a mass flow rate of 0.012 kg/s.

from 900 to 2300 W m2 K1; once vapor quality reaches the value of 0.8 the heat transfer coefficient quickly decreases because of the dry-out phenomenon. As the inlet temperature decreases, the convective heat transfer coefficient profile rises, due to the variation of

the refrigerant properties; in any case the dependence of the heat transfer coefficient on the temperature is limited. In the second test, the inlet temperature is fixed and the refrigerant mass flow rate has been varied; five different profiles are

126

M. De Carli et al. / Energy Conversion and Management 95 (2015) 120–130 Table 1 General input data for the calculations in the simplified model for the GHE.

convecve heat transfer coefficient G = 88.5

G = 132.7

G = 176.9

G = 221.2

G = 265.4

4500 4000

α (W/(m2 K))

3500 3000 2500 2000

Refrigerant

R 410A

Maximum diameter (m) Borehole diameter (m) Pipe diameters (m) Undisturbed ground temperature (°C) Ground thermal conductivity (W m1 K1) Filling thermal conductivity (W m1 K1)

10 0.15 Outer: 0.012 inner: 0.010 13 2 2

1500 1000 500 0 0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

quality Fig. 6. Profiles of the convective heat transfer coefficient for R410A and for variable specific mass flow rate (kg/(m2 s)).

reported (Fig. 6) for an inlet pressure of 799 kPa and an inlet temperature of 0 °C. The profiles of the convective heat transfer coefficient rise until a value of 0.8 approximately (where the maximum value ranges from 1500 to 4000 W m2 K1) with a slope which is different in each case; in the last part of the evaporating process the convective heat transfer coefficient quickly decreases. It has to be underlined that in this case the length of the pipe in the ground has been changed in order to allow the complete evaporation of the refrigerant fluid. The difference between two adjacent profiles is caused both by the contribution of the forced two-phase convection and by the temperature drop. As the mass flow rate grows, considering a fixed pipe diameter, Reynolds number increases, as well as the forced two-phase convection: this fact causes higher pressure losses and a wider temperature drop, as shown in Fig. 7. If the refrigerant used is R134a considerations are similar to the ones carried out for R410A case.

4.2. Simplified model of the DX-GHE In the second set of tests, the calculations of the evaporating process have been coupled with the simplified model of the ground heat exchanger in steady state conditions (Fig. 1, Eqs. (1) and (2)). There are many parameters that affect the pipe behavior, such as geometry of the borehole, drill length, soil and filling properties, refrigerant mass flow rate and its inlet temperature. Geometry and characteristics of the borehole and the soil are fixed (Table 1), while the remaining parameters vary one by one. Table 2 shows

Saturation pressure G = 88.5

G = 132.7

G = 176.9

G = 221.2

G = 265.4

900

Pressure [kPa]

800 700 600 500

the outputs for different tests carried out by changing different inputs. Varying the mass flow rate from 0.008 to 0.016 kg s1 (from the first to the third column in Table 2), results show that the evaporation is complete with a high level of superheating for the minimum mass flow rate (column 1 of Table 2), while, for the other two values of mass flow rate the process ends with x < 1. The trends of the convective heat transfer coefficient and the ones of the pressure drop are in agreement with the results of the previous section. The growth of mass flow rate requires more heat-absorption rate to complete the evaporation; if the length is not sufficient the process can be incomplete (in the first three cases the length of the pipe has been set equal to 110 m), even if the growth of temperature drop helps the thermal exchange. For this reason it was decided to analyze the heat transfer by varying the inlet temperature with the inputs reported in the columns from 4 to 6 of Table 2. Due to the variation in a range of 10 °C for the inlet temperature (passing from 5 °C to 5 °C), outlet conditions are sensibly different. In the fourth case evaporation and high superheating occur, while, in the sixth one the evaporating process is incomplete. Hence performance and heat-absorption rate, considering the evaporating process, grow as the inlet temperature decreases; this is mainly due to the difference between the fluid temperature and the undisturbed ground temperature (13 °C). The fourth and the fifth cases have similar losses and they are almost double compared to the sixth case, where the refrigerant fluid leaves the GHE with low vapor quality conditions (about 0.7). High heat fluxes and good performance may be achieved by working on mass flow rate and inlet temperature, by combining their effects. Further tests have been carried out, imposing an inlet quality equal to 0.30 and other input parameters listed in columns from 7 to 9 of Table 2. In these three cases the parametric analysis has been carried out by changing the length of the pipe, trying to maintain unchanged outlet conditions. Looking at the results of Table 2, the growth of the length of the pipe may lead to complete the evaporating process, while pressure loss and temperature drop increase appreciably. The last tests comparison refers to the outlet conditions between R410A and R134a with fixed inlet data: mass flow rate 0.013 kg s1, inlet quality x = 0.3, inlet temperature 5 °C, inlet pressures 934 kPa for R410A and 350 kPa for R134a, with an overall pipe length of 100 m. At the end of the process, for R410A x = 0.77 while for R134a x = 0.93, since R134a presents higher pressure losses than R410A (Figs. 8 and 9), therefore the temperature leap available between undisturbed soil and refrigerant increases as well as the specific heat flux. 4.3. The analytical model of DX-GHE

400 300 0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Quality Fig. 7. Profile of the saturation pressure inside the ground heat exchanger for R410A and for variable specific mass flow rate (kg/(m2 s)).

The further step of the research was to build up a model which couples the thermal model of the refrigerant fluid (for evaluating the convective heat transfer coefficient, the pressure and the temperature losses) and the GHE analytical thermal model [13]. The mean thermal resistance of the refrigerant in the evaporating process is around 103 m K W1 against 101 m K W1 for the soil and

127

M. De Carli et al. / Energy Conversion and Management 95 (2015) 120–130 Table 2 GHE outputs for the sensitivity analysis based on different inlet mass flow rate, inlet conditions and borehole lengths for the evaporating process. Case 1

Mass flow rate (kg s ) Inlet temperature (°C) Inlet pressure (kPa) U-tube length (m) Outlet conditionsa Evaporation length (m) Convective heat transfer coefficient (W m2 K1) Temperature at the end of the evaporating process (°C) Temperature drop in the evaporation (°C) Pressure at the end of the evaporating process (kPa) Pressure drop during evaporation (kPa) Outlet temperature (°C) Outlet pressure (kPa) Heat absorption rate (W) Specific heat exchanged by the pipe (W m1) a

1

2

3

4

5

6

7

8

9

0.008 5 934 110 O 98 1564 3.04 1.96 879.9 54.1 10.98 874.8 1210 12.3

0.012 5 934 110 T (x = 0.85) – 2173 1.95 3.05 851.1 82.9 1.95 851.1 1352 12.3

0.016 5 934 110 T (x = 0.74) – 2669 0.67 4.33 817.6 116.4 0.67 817.6 1431 13.0

0.014 5 680 100 O 78.5 2675 12.15 7.15 529.9 150.1 9.77 494.9 2411 24.1

0.014 0 799 100 V (x  1) 100 2588 6.86 6.86 638.1 160.1 6.86 638.1 2020 20.2

0.014 5 934 100 T (x = 0.73) – 2385 2.02 2.98 853.0 81.0 2.02 853.0 1238 12.4

0.012 0 799 70 T (x = 0.82) – 2234 2.18 2.18 746.3 52.7 2.18 746.3 1303 18.62

0. 012 0 799 80 T (x = 0.9) – 2272 3.01 3.01 726.4 72.6 3.01 726.4 1516 18.95

0. 012 0 799 90 V (x 1) 90 2271 4.58 4.58 688.9 110.1 4.58 688.9 1738 19.32

O: superheated vapor; T: two phase mix; V: vapor.

Profiles of temperature, pressure and specific heat flux in the GHE for R410A

6 pressure

900

4

800

2

temperature

700

0

600 -2

500

-4

400 300

-6

specific heat flux

200 0

20

40

Temperature [ºC]

Pressure [kPa] Specific heat flux [W/m2]

1000

60

80

-8 100

length of GHE [m] Fig. 8. Profiles for the fluids R410A along the pipe: pressure and specific heat flux (main axis), temperature (secondary axis).

101 to 102 m K W1 for the filling material. Due to this consideration the analytical method [13] can be considered coupled with a simplified refrigerant model, neglecting the parameters reliance on the vapor quality for the whole process. A mean convective heat transfer coefficient can be estimated externally the procedure, by using a suitable calculation tool for the refrigerant fluid properties, e.g. REFPROP [28]. The model of the pressure losses regards only the friction effect, for the two phase mixture at mean quality (Eq. (12)), and the motion variation between the inlet and the outlet

Profiles of temperature, pressure and specific heat flux in the GHE for R134a

900

6 4

temperature

800

2

700

0

600

specific heat flux

500

-2 -4

400 pressure

300 200 0

20

40

60

80

Temperature [ ºC]

Pressure [kPa] Specific heat flux [W/m2]

1000

-6 -8 100

length of GHE [m] Fig. 9. Profiles for the fluids R134a along the pipe: pressure and specific heat flux (main axis), temperature (secondary axis).

conditions (Eq. (17)). Gravity effect involves a gain of pressure for the down-stream and an effective loss for the up-stream, therefore, from the inlet to the outlet, the pressure variation due to gravity can be neglected. The pressure losses model depends on the inlet/outlet conditions and on the mean vapor quality. Since the mass flow rate is fixed (Eq. (18)), as well as the inlet quality and the numbers of drills, the analytical method can be applied supposing constant the refrigerant temperature and neglecting the ground penalty temperature in a first stage. The length obtained in this way can be used to find the pressure losses and then the outlet temperature at the GHE. After few iterations the length and the outlet temperatures are determined. Once sized the overall length of the GHE, the ground penalty temperature can be calculated, thus correcting the size of the GHE. The proposed simplified model has been run, based on a typical profile of the heating energy demand of a residential building (Table 3). Table 4 summarizes the results of the simplified method. As can be seen, if the ground temperature penalty is considered, the length grows of about 5%. 4.4. Overall model of the DX-GSHP cycle In this section the model includes the calculation of the GHE based on the analytical ASHRAE method as well as the heat pump behavior. Results will refer to R410A, but they are similar to R134a. Also in this case the typical heating energy need of the residential building of Table 3 has been considered. The heating capacity of the heat pump is 8.8 kW with a nominal COP of 3.8. The code needs the number of boreholes, inlet vapor quality, temperature and pressure. As a result of the simulation, the sizing of the ground heat exchangers is carried out by searching the complete evaporation

Table 3 Typical mean monthly heating energy needs for a residential building in Venice (kW h). January February March April May June July August September October November December

1971 1481 932 157 0 0 0 0 0 159 1059 1822

128

M. De Carli et al. / Energy Conversion and Management 95 (2015) 120–130

Table 4 Most relevant results of the DX-GSHP sizing with the analytical model. Mass flow rate for each GHE (kg s1) Inlet quality Undisturbed ground temperature (°C) Inlet temperature (°C) Inlet pressure (kPa) Pressure drop (kPa) Temperature drop (°C) Meters depth required for each GHE without penalty temperature (m) Total meters depth required for the whole system without penalty temperature (m) Mean convective heat transfer coefficient (W m2 K1) Meters depth required for each GHE with penalty temperature (m) Total meters depth required for the whole system with penalty temperature (m)

0.015 0.3 13.1 5 934 127.2 6.55 77 231 2850 81 243

outlet condition varying the total length of the U-tube. As for the boreholes field, three boreholes coupled in parallel and arranged along a line spaced 7 m apart are considered. Table 5 summarizes the main results of the sizing calculation; the convective heat transfer coefficient, the pressure and the temperature profiles inside the pipe are reported in Fig. 10. As the program allows the calculation of the refrigerant thermodynamic cycle, using the current energy needs it is possible to estimate the effective COP. In this way it is possible to compare the design COP of the heat pump (that is 3.8) with the value resulting from simulations (that is 3.7). The results of the calculations are listed in Table 6, together with outputs of a further simulation which has been carried out by

Mass flow rate for each GHE (kg s1) Inlet quality Undisturbed ground temperature (°C) Inlet temperature (°C) Inlet pressure (kPa) Meters depth required for each GHE (m) Total meters depth required for the whole system (m) Mean convective heat transfer coefficient (W m2 K1) Pressure drop (kPa) Temperature drop (°C) Heat absorption rate (W) Performance of the borehole (W/m) Condensing temperature (°C) Compressor electric power (W) Condenser thermal power post-sizing (W) Resulting COP

0.015 0.3 13.1 5 934 65.5 196.5 2586 238.1 9.1 6570 33.2 40 2399 8970 3.7

Pressure [kPa] Heat exchange coefficient [W/(m2 K)]

Table 5 Most relevant results of the DX-GSHP sizing.

6

3500 heat exchange coefficient

3000

5

3

temperature

2

2000

1 1500

0 pressure

-1

1000

Temperature [ ºC]

4 2500

-2 500

-3 -4

0 0

20

40

60

80

100

120

distance from the inlet of the U-tube Fig. 10. Convective heat transfer coefficient, pressure and temperature profiles in the considered evaporating process.

Table 6 Sizing and post correction of DX-GSHP systems with different inlet temperature at the GHEs. Mass flow rate for each GHE (kg s1) Inlet quality Undisturbed ground temperature (°C) Inlet temperature (°C) Inlet pressure (kPa) Meters depth required for each GHE (m) Total meters depth required for the whole system (m) Mean convective heat transfer coefficient (W m2 K1) Pressure drop (kPa) Temperature drop (°C) Heat absorption rate (W) Performance on borehole meter (W m1) Condensing temperature (°C) Compressor electric power (W) Condenser thermal power post-sizing (W) Resulting COP

0.015

0.015 0.3 13.14

5 934 65 195 2598 221. 8.49 6444 33.0 40 2372 8816 3.72

0 799 45.5 136.5 2595 175 7.43 6519 47.8 40 2435 8953 3.67

reducing the starting evaporating temperature, maintaining all the other parameters unchanged. The length sizing decreases and the borehole performance increases as the inlet temperature is reduced, but the COP slightly decreases; as a matter of fact, the increase of the length causes higher pressure and temperature drops, thus penalizing the compressor work. There are many parameters that can be varied to compare the different results like the mass flow rate, the borehole length, the quality, the drills number and the inlet temperature; the last one is the most suitable relating the approach used and it has a large influence on the GHEs sizing length, as well as on the heat pump COP. The results of the overall model of the DX-GSHP can be compared with the analytical method without considering the heat pump cycle. A comparison can be carried out considering Tables 4 and 5, which show the results based on the same input and operating conditions of the heat pump. With the simplified analytical method temperature and pressure losses are lower and, due to this result, the required length is oversized by about 24%, achieving a conservative result. 4.5. Comparison between DX-GSHP and SL-GSHP Finally, a comparison between a direct expansion system (DXGSHP) and a traditional one (SL-GSHP) has been carried out. For this purpose, the same energy needs (Table 3) and design conditions have been used. The heat pump of the SL-GSHP has the same thermodynamic cycle of the DX one. The nominal heating capacity of the heat pump is 8.8 kW in both cases. As for the ground heat exchangers the lay-out for both the systems consists of three boreholes connected in parallel. As for the GHEs, grouting diameter is 0.15 m for the DX-GSHP, while it is 0.13 m in the case of the SLGSHP. The pipe of the DX-GSHP has 10 mm as inner diameter and 12 mm as outer diameter. For the SL-GSHP a single U-tube has been considered with internal and external diameters of 26.2 mm and 32 mm, respectively. Table 7 summarizes the results of the comparison between the two systems. In Case 1 the comparison has been carried out in order to maintain the same COP. The length required by the traditional system is appreciably higher than that of the direct expansion heat exchanger (265 m versus 136 m). As for the costs, it has to be underlined that, even if the depth is lower in the DX-GSHP copper pipes are more expensive than plastic pipes, therefore the installation cost analysis has to be carried out carefully. A further test (Case 2) has been carried out trying to obtain the same mean temperature of the fluid inside the GHEs for both SL and DX heat pumps. A temperature of R410A around 3.5 °C and 4 °C have been considered for the DX and the SL systems

129

M. De Carli et al. / Energy Conversion and Management 95 (2015) 120–130 Table 7 Most relevant data for the comparison between the DX-GSHP and the SL-GSHP in the different cases. Case 1

Design COP Inlet temperature at the GHE (°C) Mean refrigerant temperature (°C) Boreholes number Total meters required for the drills (m) Pressure losses in the GHEs system (kPa) Compressor electric power (W) Power of the auxiliary system (W)

Case 2

Case 3

Case 4

DX

SL

DX

SL

DX

SL

DX

SL

3.66 0 3.4 3 136 175 2435 –

3.52 2 4 3 265 68 2411 86

3.75 9 3.5 3 288 326 2386 –

3.52 2 4 3 265 68 2411 86

3.66 0 3.4 3 136 175 2435 –

3.95 2 4 3 292 99 2100 126

4.20 5 1.4 3 216 190 2115 –

3.95 2 4 3 292 99 2100 126

respectively. In this case the DX length sizing leads to a GHE longer than the SL one, about 23 m. The difference is appreciably, ideally this value may be lower considering the input data and the working conditions of the heat pump. Due to the different COP between the systems, in the direct expansion heat exchanger more thermal load has to be absorbed by the ground. Furthermore, the operations are different: for sizing the SL-CHP the mean typical parameters of the process are used, while, for the DX-CHP they are variable along the whole U-tube. This fact may explain the difference of the length calculated for the DX and the SL systems. On the reported tests, the compressor has the same feeding power for both the systems. Traditional systems have, in the real case, less compressor power than the supposed one, due to the pressure losses realized by the evaporating fluid within the water/refrigerant heat exchanger that are usually around 10 kPa. This value involves a decrease of the saturation temperature of about 1.5 °C to 2 °C for the fluid R410A at the usual working pressure conditions. It has to be underlined that, if the working pressure decreases, the losses cause an increase in the temperature drop, due to the divergence of the isobars. This may explain the above reported values. In Case 3 the comparison is carried out adjusting the compressor power and considering also the pressure losses of the water in the evaporator of the heat pump of the SL system. In this case the total pressure losses (in the evaporator of the heat pump as well as in the horizontal loop and in the GHEs) of the auxiliary pump are considered. As for the compressor, the compression process is supposed ideal and R410A mass flow rate is calculated at the evaporator using the first principle of thermodynamic and considering only the evaporating process (Dx = 0.8 and a mass flow rate of 0.0385 kg s1). The compressor power decreases about 400 W and the auxiliary power grows. As a result, the whole system presents a higher COP compared to Case 2 for the SL-GSHP. Another important result is the total depth required for the drills: DX system allows to size the GHEs with a lower length than the SL one. In the last comparison (Case 4) the inlet parameters change: inlet values are x = 0.2, t = 5 °C and a mass flow rate 0.0138 kg s1 for each drill. The DX-GSHP presents a higher COP and less required length of the GHEs compared to the SL- system. 5. Conclusions The present study shows an analysis on the DX-GSHP, which considers the thermodynamic process of the refrigerant in the evaporation through the GHE. For this purpose a mathematical model has been developed. The program allows to model the two-phase flow (calculation of the quality, convective heat transfer coefficient, temperature and pressure losses) coupled with the heat transfer through the ground, based on an analytical method. The vapor quality is considered as an independent variable, which allows an easy link between the two-phase outflow model and the calculation of the heat transfer through the ground.

The thermal process, and hence the convective heat transfer profile, varies appreciably with respect to the inlet quality, as well as to the inlet temperature and to the specific mass flow rate. These three parameters have a mutual influence on the heat transfer coefficient. The analysis on the refrigerant outflow allows to calculate the temperature drop that has a great influence on the heat exchange with the ground. Once checked the combination of the model of the heat transfer with the ground with the model of the evaporating process, the heat pump has been included, thus allowing the simulation of the entire DX-GSHP system. The complete model has been used for a comparison with a traditional water based system (SL-GSHP), showing that, using the same heating energy demands and the same design heating capacity of the heat pump, the COP of the DX-GSHP is higher (in terms of COP in a range between 0.15 and 0.25) and the length of the heat exchangers is lower than that obtained with the SL-GSHP. If the comparison is carried out considering the same mean fluid evaporating temperature inside the GHEs, the COP of the DX system is still slightly higher than the SL system, but the overall length of the GHE for the DX system is greater than the one of the water based system. The model carried out in this work is relatively flexible in terms of input conditions, such as heat pump and soil characteristics, borehole dimensions and patterns, as well as initial evaporating conditions inside the GHE. Further experimental investigation should deal with the seasonal performance of a DX-GSHP and a SL-GSHP in terms of costs and Total Equivalent Warming Impact (TEWI), considering both the energy consumed by the systems and the overall impact in terms of equivalent CO2 emissions. References [1] Shi L, Chew MYL. A review on sustainable design of renewable energy systems. Renew Sustain Energy Rev 2012;16(1):192–207. [2] Bayer P, Saner D, Bolay S, Rybach L, Blum P. Greenhouse gas emission savings of ground source heat pump systems in Europe: a review. Renew Sustain Energy Rev 2012;16:1256–67. [3] Li M, Lai ACK. Thermodynamic optimization of ground heat exchangers with single U-tube by entropy generation minimization method. Energy Convers Manage 2013;65:133–9. [4] Vakiloroaya V, Samali B, Fakhar A, Pishghadam K. A review of different strategies for HVAC energy saving. Energy Convers Manage 2014;77:738–54. [5] Wang X, Ma C, Lu Y. An experimental study of a direct expansion groundcoupled heat pump system in heating mode. Int J Energy Res 2009;33:1367–83. [6] Wang H, Zhao Q, Wu J, Yang B, Chen Z. Experimental investigation on the operation performance of a direct expansion ground source heat pump system for space heating. Energy Build 2013;61:349–55. [7] Fannou JL, Rousseau C, Lamarche L, Stanislaw K. Experimental analysis of a direct expansion geothermal heat pump in heating mode. Energy Build 2014;75:290–300. [8] Beauchamp B, Lamarche L, Kajl S. A dynamic model of vertical direct expansion ground heat exchanger. In: Proceedings of international conference on renewable energies and power quality 2008, Santander, Spain; 2008. [9] Austin BT, Sumathy K. Parametric study on the performance of a directexpansion geothermal heat pump using carbon dioxide. Appl Therm Eng 2011;31:3774–82.

130

M. De Carli et al. / Energy Conversion and Management 95 (2015) 120–130

[10] Eslami-Nejad P, Ouzzane M, Aidoun Z. Modeling of a two-phase CO2-filled vertical borehole for geothermal heat pump applications. Appl Energy 2014;114:611–20. [11] Ingersoll LR, Zobel OJ, Ingersoll AC. Heat conduction: with engineering and geological applications. New York: McGraw-Hill Book Co.; 1954. [12] Carslaw HS, Jaeger JC. Conduction of heat in solids. Oxford: Claremore Press; 1959. [13] Kavanaugh SP, Rafferty K. Ground source heat pumps - design of geothermal systems for commercial and institutional buildings. ASHRAE Applications Handbook; 1997. [14] Eskilson P. Thermal analysis of heat extraction boreholes. Doctoral thesis, Department of Mathematical Physics and Building Technology, University of Lund, Lund, Sweden; 1987. [15] Hellström G, Sanner B. Earth energy designer: software for dimensioning of deep boreholes for heat extraction. Department of Mathematical Physics, University of Lund, Lund, Sweden; 1994. [16] Kurevija T, Vulin D, Krapec V. Effect of borehole array geometry and thermal interferences on geothermal heat pump system. Energy Convers Manage 2012;60:134–42. [17] Zarrella A, Scarpa M, De Carli M. Short time-step performances of coaxial and double U-tube borehole heat exchangers: modeling and measurements. HVAC&R Res. 2011;17(6):959–76. [18] Zarrella A, Scarpa M, De Carli M. Short time step analysis of vertical groundcoupled heat exchangers: the approach of CaRM. Renewable Energy 2011;36:2357–67.

[19] Paul ND. The effect of grout thermal conductivity on vertical geothermal heat exchanger design and performance, Master’s thesis, South Dakota State University, US; 1996. [20] Remund CP. Borehole thermal resistance: laboratory and field studies. ASHRAE Trans 1999;105(1). [21] Watson KM. Thermodynamics of the liquid state – generalized prediction of properties. Ind Eng Chem 1943;35:398–406. [22] Collier JC. Boiling and evaporation. In: Heat exchanger design handbook. Hemisphere; 1983. [23] Chen JC. Correlation for boiling heat transfer to saturated liquids in convective flow. Ind Eng Chem Process Des Dev 1966;5:322. [24] Hewitt GF, Shires GL, Bott TR. Process heat transfer. Boca Raton: CRC Press Begell House; 1994. [25] Blasius H. Das Aehnlichkeitsgesetz bei Reibungsvorgängen in Flüssigkeiten. Vereines deutscher Ingenieure. Mitteilungen über Forschungsarbeiten auf dem Gebiete des Ingenieurwesens 1913;131. [26] McAdams WH, Woods WK, Heroman Jr LC. Vaporization inside horizontal tubes-II-benzene–oil mixtures. ASME Trans 1942;64:193. [27] Capozza A, De Carli M, Zarrella A. Design of borehole heat exchangers for ground-source heat pump: a literature review, methodology comparison and analysis on the penalty temperature. Energy Build 2012;55(12): 369–79. [28] National Institute of Standards and Technology-NIST. Reference fluid thermodynamic and transport properties database (REFPROP): Version 7.