Influence of an internal heat exchanger on the operation of a CO2 direct expansion ground source heat pump

Energy & Buildings 202 (2019) 109343

Contents lists available at ScienceDirect

Energy & Buildings journal homepage: www.elsevier.com/locate/enbuild

Influence of an internal heat exchanger on the operation of a CO2 direct expansion ground source heat pump A. Nguyen a,∗, P. Eslami-Nejad a, M. Badache a, A. Bastani a CanmetENERGY, Naturel Resources Canada, 1615 Lionel-Boulet, P.O. Box 4800, Varennes, Québec J3X 1S6, Canada

a r t i c l e

i n f o

Article history: Received 1 May 2019 Revised 22 July 2019 Accepted 29 July 2019 Available online 29 July 2019 Keywords: DX-GSHP Borehole heat exchanger Internal heat exchanger Ground surface temperature CO2 Transient simulation

a b s t r a c t In this work, the influence of an internal heat exchanger (IHE) on the operation of a CO2 direct expansion ground source heat pump (DX-GSHP) was assessed through transient energy simulations. A previous theoretical model was extended to account for an IHE, and was used to evaluate the performance of DXGSHP systems under various controls, ground conditions and building loads. A showcase of the benefit of integrating an IHE to the DX-GSHP was carried out, through 4-month heating energy simulations period of a small residential building located in Canada, under 4 different cases. Results show that adding an IHE to the system effectively counters the effect of thermal short-circuiting in the borehole as well as ground surface temperature variations, and thereby increases its efficiency.

1. Introduction It is well established that ground source heat pump (GSHP) systems are an environmentally friendly technology that have the potential to greatly reduce the overall energy consumption of buildings. Over the years, they have been adopted widely around the world due to their high efficiency, low maintenance cost and low green house gas emissions. These benefits are achieved through the borehole heat exchanger (BHE), which exchanges heat with the surrounding ground and serves as a heat source/sink for the system. Among all the existing types and configurations of BHEs, the secondary-loop vertical U-tube is by far the most common type of BHE [26]. For decades, numerous studies have been conducted on a variety of subjects (e.g. model developments, parametric analysis, simulations, control sequences, design methods, new types of BHE) to foster the use of GSHP systems worldwide. Despite these efforts, GSHP systems are still faced with major market barriers: high initial cost and long pay back periods [14]. To tackle these barriers, developments focusing on increasing the overall performance of the system while reducing the cost associated with the underground components are mandatory. Recently, the less common direct-expansion ground source heat pump (DX-GSHP) has gained some research attention due to the expected emergence of CO2 as refrigerant [29]. Naturally, CO2 is ∗

Corresponding author. E-mail address: [email protected] (A. Nguyen).

https://doi.org/10.1016/j.enbuild.2019.109343 0378-7788/Crown Copyright © 2019 Published by Elsevier B.V. All rights reserved.

Crown Copyright © 2019 Published by Elsevier B.V. All rights reserved.

convenient for DX-GSHP applications since it is non-toxic, nonflammable, non-corrosive, accessible and inexpensive. Also, its immiscibility with mineral oil make CO2 a great choice to mitigate the oil return problem [18]. Although CO2 DX-GSHP systems have the potential to be more economical then conventional ones [12], too little researches have been conducted up to now as it is a relatively new subject of interest. In fact, only a handful of models have been developed for CO2 DX-GSHPs, all with different degrees of complexity, from basic steady-state models to more complex transient models. Austin and Sumathy [1] proposed a steady state model of a CO2 DXGSHP in a transcritical cycle where they investigated the relationship between gas cooler and evaporator sizing for horizontal loops. Eslami-Nejad et al. [9] developed a quasi-transient model for a CO2 DX-GSHP in a transcritical cycle for heating applications with DX-BHEs. They showed that, due to changes in the system working conditions, the total borehole length is not linearly correlated with performance parameters such as heat pump energy consumption and coefficient of performance, and pressure drop in the borehole. Ghazizade-Ahsaee and Ameri [11] did a similar work to Austin and Sumathy [1] and studied the steady-state performance of a transcritical direct-expansion CO2 GSHP with horizontal loops for space-heating and DHW applications. Variable speed compression in CO2 DX-GSHP was studied numerically by [20]. They extended the work of [19] and developed a transient coupled model of a variable speed CO2 DX-GSHP using a basic compression cycle for space heating and cooling. They showed through annual energy simulations that, in heating mode, the BHE cannot efficiently

2

A. Nguyen, P. Eslami-Nejad and M. Badache et al. / Energy & Buildings 202 (2019) 109343

Nomenclature c C COP d D g f fr G H k L m˙ M n P Q qg r R t td T UA v Vswept W z

specific heat capacity (J kg−1 ◦ C−1 ) thermal capacity (J ◦ C−1 ) coefficient of performance (-) buried depth (m) shank spacing (m) gravitational acceleration (m s−2 ) friction factor (-) frequency ratio (-) mass flux (kg s−1 m−2 ) enthalpy (J kg−1 ) thermal conductivity (W m−1 ◦ C−1 ) BHE length (m) mass flow rate (kg s−1 ) mass (kg) number nodes (-) pressure (Pa) heat transfer rate (W) geothermal gradient (◦ C m−1 ) radius (m) thermal resistance (m ◦ C W−1 ) time (s) time phase difference (s) temperature (◦ C) overall heat transfer coefficient (W ◦ C−1 ) volume (m) swept volume rate (m3 s−1 ) compressor work (W) depth (m)

Greek symbols α thermal diffusivity (m2 s−1 ) ηis isentropic efficiency (-) ηv volumetric efficiency (-) μ viscosity (Pa s) ρ density (kg m−3 ) ω angular frequency (s−1 ) Subscripts 1 suction of compressor 2 discharge of compressor 3 outlet of condensor/gas cooler 4 outlet of internal heat exchanger (high pressure) 5 inlet of evaporator 6 inlet of internal heat exchanger (low pressure) a amplitude AUX auxiliary b bulk or borehole f fluid g grout HP heat pump i time step index j node index k neighbouring node index m mean p pipe s soil w wall z vertical direction Acronyms BHE borehole heat exchanger DX direct-expansion IHE internal heat exchanger

LSHE GSHP GSTV TRCM

load-side heat exchanger ground source heat pump ground surface temperature variation thermal resistance and capacity model

provide the required superheat at its outlet due to the existence of a strong thermal short-circuiting around the ground surface when the system is operating at low mass flow. It is worth noting that in all models cited above, none accounted for the effect of seasonal ground surface temperature variation (GSTV). In reality, getting the required amount of superheat at the outlet of the BHE, in heating mode, can be even more difficult as the ground temperature near the surface is significantly colder. While it is recognized that GSTV affects the performance of conventional borehole heat exchanger fields [5,7,13,16,23,27], it is unclear how much it can impact the performance of direct expansion BHEs, especially short ones (i.e. below 45 m of depth). A possible solution to counter GSTV and thermal shortcircuiting is to maintain two-phase flow in the BHE by the addition of an internal heat exchanger (IHE) to the cycle, as suggested by [20]. This way, the required superheat would be provided by the IHE, rather than by the BHE. This was analyzed by [10] for horizontal CO2 DX-GSHPs in steady-state operation. Their results showed that the IHE procured only a slight benefit to the system. However, the ground heat exchanger in their model considered a uniform heat flux and a fixed temperature at the radius 1 m, thereby not accounting for thermal short-circuiting and transient ground behavior. The works presented hereinafter are a direct extension of previous works of [20] and aim to present a CO2 DX-GSHP model that accounts for the complete refrigerant cycle using an IHE and assess the influence of GSTV in transient energy simulations. The aims of this work are twofold: 1) present the methodology for modeling the DX-GSHP and 2) showcase the benefit of integrating an IHE to the DX-GSHP model through simulations. 2. Modeling approach The transient modeling of a DX-GSHP requires the evaluation of each point of the refrigerant cycle in terms of pressure and enthalpy at every time step. For basic parameter evaluation of the fluid, the open source library CoolProp [4] is used. The model is programmed in a compact script and integrated with a commercial ODE solver provided by MATLAB [28]. This section is subdivided, for clarity, into three parts: 1) mathematical relationships governing the different heat pump components, 2) ground surface temperature variations and 3) simulation strategy. 2.1. Heat pump components The DX-GSHP system under this study is composed of four main components; a compressor, an expansion valve and two heat exchangers: one on the sink side to provide heating to the building (load-side heat exchanger) and one on the source side (here a vertical DX-BHE). An additional component (internal heat exchanger) is also used to improve the performance of the heat pump. The schematic of the DX-GSHP system and the corresponding vaporcompression cycle are presented in Fig. 1(a) and (b). Each component of the heat pump and the governing equations are briefly described the following subsections. 2.1.1. Compressor In order to model the compressor, volumetric (ηv ) and isentropic (ηis ) efficiency functions are required. In single speed compressors, efficiency is a function of compression ratio of the heat

A. Nguyen, P. Eslami-Nejad and M. Badache et al. / Energy & Buildings 202 (2019) 109343

3

Fig. 1. a) Schematic of the DX-GSHP in heating mode and b) Illustration of the complete transcritical cycle.

Fig. 2. Illustration of the a) isentropic efficiency and b) volumetric efficiency contour of the compressor model in this work.

pump (P2 /P1 ). However, in variable speed compressors, frequency alters the efficiency. Therefore, in this study compressor efficiencies are calculated as a function of the compression ratio, as well as the frequency ratio (fr ). This is done by normalizing the efficiency terms for a speed of 60 Hz, which is where the compressor operates with the maximum mechanical and motor efficiency. The relations are given by:

η = 0.91 + 0.18 fr − 0.09 fr .2 η60 ηis,60 = 1.003 − 0.121

P2 P1

ηv,60 = 1.079 − 0.1053

P2 P1

(1)

H2,is − H1

ηis

m˙ = ηv ρ1Vswept fr

(5)

where ρ 1 is the density at the suction of the compressor and Vswept is the swept volume rate at 60 Hz. In this work, the discharge pressure is controlled using the correlation developed by [25], given by the following expression:

P2 = 132.3 − 8.4T3 + 0.3T32 − 27.7 × 10−4 T33

(6) ◦ C,

(2) (3)

These relationships were obtained by best fitting the manufacturer data of a carbon dioxide compressor. For simplicity and to avoid modeling off/on cycles based on the indoor temperature control and temperature deadbands, it is assumed that the speed can vary from 0 to 66 Hz. Both the contour of ηis and ηv are depicted in Fig. 2. The fluid enthalpy after compression is given by:

H2 = H1 +

nally, the mass flow in the heat pump is calculated through:

(4)

where H1 is the enthalpy of the fluid at the suction of the compressor and H2,is is the isentropic enthalpy after compression. Fi-

where P2 and T3 are expressed in bar and respectively. Any other existing correlations can be implemented in the model; however, this correlation was chosen for its simplicity and low temperature applicability. According to [25], it is valid for T3 ranging from 25 to 45 ◦ C. 2.1.2. Load side heat exchanger In practice, fin-tube cross-flow heat exchangers are used for heat transfer between refrigerant and air. However, here, to simplify the modeling process a counter-flow tube-in-tube heat exchanger is used in which CO2 flowing in the inner-tube and air flowing through the annulus. Although a tube-in-tube heat exchanger might have different performance compared to fin-tube cross-flow heat exchanger, it provides a similar representation of CO2 temperature evolution inside the heat exchanger. The load side heat exchanger (LSHE) of the heat pump is modeled numerically

4

A. Nguyen, P. Eslami-Nejad and M. Badache et al. / Energy & Buildings 202 (2019) 109343

using a constant overall heat transfer coefficient (UA). However, the heat transfer equation is discretized over n nodes as shown by:

Qi =

UA · (T23,i−1 − Tair,i−1 ) n

T23 = f (P23 , H23 )

(7) (8)

where Q is the heat exchange rate, T is the temperature and i is the node index. The enthalpy change of CO2 along the tube is also calculated using the heat exchange rate of each node (Qi ) assuming that no heat is transferred to the environment. The energy balances for each node for both CO2 and air are expressed as:

H23,i

Qi = H23,i−1 − m˙ 23

Tair,i = Tair,i−1 −

(9)

Qi m˙ air c p,air

(10)

where P is the pressure, H is the enthalpy, m˙ is the mass flow, cp is the specific heat and i is the node index. The total heat exchange rate is thus:

Qtot = m˙ 23 (H23,1 − H23,n ) = m˙ air c p,air (Tair,n − Tair,1 )

2.1.3. Expansion device In the heating mode (the focus of this work) the expansion valve drops the LSHE pressure (high pressure side of the cycle) to BHE pressure (low pressure side). The following equation governs the expansion process based on an isenthalpic expansion:

(12)

The evaporating temperature is calculated so that the amount of superheat at T1 (5 ◦ C) is satisfied. The lowest pressure allowed for this study is set to 2 MPa (−19.5 ◦ C). If this pressure is reached, the heat pump will cease and the auxiliary heater will provide heating to the building. 2.1.4. DX borehole heat exchanger The DX-BHE model accounts for both the BHE and the ground. The quasi-3D thermal resistance and capacity approach, known as the thermal resistance and capacity model (TRCM) [3,19,21,22,24,30] is used and represented schematically in Fig. 3. The model accounts for the enthalpy change of the fluid through heat transfer by advection and conduction in the fluid (CO2 ) and in the adjacent grout and ground respectively. For the pipe, grout and ground, the time derivative of the temperature T for a given control volume is given by: nj

Cj

 Tk − T j dT j = dt Rk k=1

∀ j = 1...n

Mj

Tpi, j − T f, j dH j H = + dt Rf Rz

(13)

∀ j = 1 . . . nz

(14)

where Tf is the fluid’s bulk temperature, Tpi is the temperature of the inner pipe wall, M is the mass associated to each control volume of the fluid, Rf is the fluid thermal resistance and Rz is the vertical advective resistance connecting the fluid nodes. In Eq. (14), the first term on the right hand side represents the enthalpy change due to the effect of heat transfer and the second term represents enthalpy change due to advection. The model also accounts for changes in fluid pressure as it is coupled to its temperature. The fluid pressure is calculated using the following steady state momentum equation:

(11)

where the index 1 denotes the inlet of the CO2 and n denotes the inlet of air. Here, the inlet air temperature is kept constant at 22 ◦ C throughout the year while the supply air temperature is set to 36 ◦ C. For simplicity, it is assumed that CO2 pressure drop is insignificant. Since the building load at each time step is known, it is convenient to solve the model in regards to the load. Hence, the solution starts with input value of Tair,1 , TCO2,1 , m˙ air (estimated based on the building load) and a guess value of m˙ CO2 . Using the above equations, the subsequent node temperatures are sequentially calculated. Then, m˙ CO2 is iteratively updated (using the secant method for instance) until the calculated inlet air temperature (e.g. Tair,n ) matched the input inlet air temperature within the acceptable tolerance. The best result in terms of precision and computation was obtained with n = 10 0 0, which assured required precision for the whole simulation period.

H4 = H5

where C is the thermal capacity, T is the temperature, R is the thermal resistance between the connected nodes, n is the total number of nodes in the network, j is the node index, nj is the number of neighbouring nodes to node j and k is the index of the neighbouring node. Note that vertical heat transfer in the ground is accounted for. The temporal derivative of the fluid’s enthalpy H for a given control volume is governed by the rate of internal energy change:

P (t, j ) = P0 (t ) −

j 

Pk (t )

∀ j = 2 . . . nz

(15)

k=1

where P0 (t) is the pressure at the inlet of the BHE and Pk (t) is the pressure change across the control volume at a given time. The latter is a function of several factors including: frictional pressure drop (Pf ), acceleration pressure drop (Pa ) and gravitational pressure drop (Pg ). Each term can be calculated using the following above equations:

P = Pf + Pa + Pg

(16)

G2 dz 4r pi ρ¯

(17)

P f = f

 1  1 Pa = G2 · − ρout ρin

(18)

Pg = ±gρ¯ dz

(19)

where f is the friction factor, G is the fluid’s mass flux, g is the gravitational acceleration, ρ¯ is the mean density of the fluid in the control volume, and ρ in and ρ out are the fluid’s density entering and exiting the control volume. All the required details for Eqs. (13)–(19) are given in [19]. 2.1.5. Internal heat exchanger The same approach as the one employed for the LSHE is used to model the internal heat exchanger (IHE), except in this case CO2 is flowing in both tube and the annulus of the heat exchanger. In short, CO2 at the condenser/gas cooler outlet is flowing inside the tube of the IHE (34 index) and CO2 at the evaporator outlet is entering into the annulus of the IHE from the opposite side (16 index). Similar to LSHE energy balance equations are written as follows:

Qi =

UA · (T34,i−1 − T16,i−1 ) n

T34,i = f (P34,i , H34,i ) H34,i = H34,i−1 −

Qi m˙ 34

T16,i = f (P16,i , H16,i ) H16,i = H16,i−1 −

Qi m˙ 16

Qtot = m˙ 34 (H34,1 − H34,n ) = m˙ 16 (H16,1 − H16,n )

(20) (21) (22) (23) (24) (25)

A. Nguyen, P. Eslami-Nejad and M. Badache et al. / Energy & Buildings 202 (2019) 109343

5

Fig. 3. Illustration of the TRCM network with 3 layers.

Recall that the evaporating pressure is adjusted at each time step so that the amount of superheat at T1 is exactly 5 ◦ C. 2.2. Ground surface temperature variations The mean ground temperature results from the thermal equilibrium between the ambient air temperature, surface solar radiation and the geothermal heat flux of the Earth’s crust. Near the surface, however, the ground temperature varies over the seasons and is mostly dependant on the intensity of solar radiation, which depends on multiple factors such as the distance from the Sun, solar cycle, geographic location on Earth, land morphology, weather conditions, etc. The data regarding GSTV is not always available (means limited in space and in time) and thus, various ground temperature models were developed over the years, from

simple analytical model [8,15] to more advanced numerical models [17]. To simplify the required computation of ground temperature [2] recently proposed a hybrid model where energy balance equation at the ground surface is supplemented by an empirical correlation. For simplicity, the ground temperature profile in this work is based on the well-known analytical solution proposed by [15]:

 Tg (z, t ) = Tm + qg z − Ta exp −z



ω 2α



 · cos

ω (t −td )−z



ω 2α



(26) where Tm is the yearly mean ground surface temperature, z is the depth, qg is the geothermal gradient, Ta is the amplitude of the seasonal variation, ω is the angular frequency, α is the ground

6

A. Nguyen, P. Eslami-Nejad and M. Badache et al. / Energy & Buildings 202 (2019) 109343

3. Case study During the heating season, the heat pump must control the low pressure side (i.e. at the BHE) so that the required superheat can be provided at the suction of the compressor. To assess the influence of an IHE, as a potential solution to mitigate thermal shortcircuiting in the BHE, as well as influence of GSTV on the transient performance of the DX-GSHP, comparison studies between two heat pump models are conducted: one with an IHE, and one without, for cases with and without an GSTV. The four case studies are listed in Table 1. An hourly building load profile (Fig. 5) of a 4-month heating period (November 1st to March 1st) was used to emulate a typical residential building located in Canada. The total heating demand is 5805 kWh with a peak load of 4.05 kW. Note that the load only represents a fraction of a single family building to accommodate use of a single borehole system. A summary of the parameters used for the DX-GSHP is presented in Table 2. Fig. 4. Comparison of the ground temperature profile computed with Eq. (26) and the DX BHE model (far field) at different times.

thermal diffusivity, t is the time and td is the time phasing difference corresponding to the beginning of start-up time. Note that Eq. (26) indirectly accounts for solar radiation at the ground surface through seasonal surface temperature variations. To numerically implement this into the DX-BHE model, the initial ground temperature is set to the solution of Eq. (26) calculated at the beginning of the simulation time. In addition, the boundary condition at the surface of the model (i.e. first layer of Ts ) is set as a Dirichlet type and varies over time according to Eq. (26) computed at the buried depth (d) of the BHE, which in this case is a typical value of 1.2 m. Note the amplitude of the variation still reaches around 10 ◦ C at that depth. With the aim to verify the validity of the proposed method for simulating GSTV, Eq. (26) was used to generate for four (4) reference solutions using the following parameters: Tm = 10 ◦ C; qg = 0.02 ◦ Cm−1 ; Ta = 20 ◦ C; α = 1.04·10−6 m2 s−1 ; and td = 51 days, from November 1st to February 1st. Fig. 4 presents a visual appreciation of the model’s ground temperature profile at the far field and the solution of Eq. (26). As shown, results indicate very good agreement between the two models for all 4 times with an average mean absolute difference (MAD) of 0.02 ◦ C. This validation shows that the model proposed in this work can be used to emulate the ground temperature profile in response to GSTV.

2.3. Simulation strategy To begin a simulation, the DX-GSHP specifications (BHE, LSHE, IHE and compressor) and the operating parameters (return and supply air temperature set points, pressure control and amount of superheat) are input. Next, one must define the initial conditions of the model in terms of pressure, enthalpy and temperature (i.e. P(t0 ) and H(t0 ) corresponding to the initial ground temperature T(t0 )). At each time step, the ODE solver outputs H (ti+1 ) and T (ti+1 ) by solving Eqs. (13)–(14) while a sequential process is used to solve Eq. (15). The solution found ensures that the load provided by the heat pump matches the heating load, while respecting the heat pump controls, at the time step. This is done by varying the frequency of the compressor (up to the maximum allowed frequency). If the building load exceeds the capacity of the heat pump at a given time, the remaining load is provided by an electrical auxiliary system.

Table 1 Case description. Parameter

Case #1

Case #2

Case #3

Case #4

Unit

Tm qz Ta IHE

10 0.02 20 Yes

10 0.02 0 Yes

10 0.02 20 No

10 0.02 0 No

◦

C C m−1 ◦ C ◦

Table 2 Parameters used for the DX-GSHP model. Component

Parameter

Value

Unit

BHE

rpi rpo rb rs D d L kp kg ks

ρ p cp ρ g cg ρ s cs η

0.0032 0.0040 0.039 100 0.023 1.2 45 385 1.5 2 0.35 2.25 2.4 Eq. (1)

m m m m m m m W m−1 ◦ C−1 W m−1 ◦ C−1 W m−1 ◦ C−1 MJ m−3 ◦ C−1 MJ m−3 ◦ C−1 MJ m−3 ◦ C−1 -

Compressor

fr Vswept P2 UA

0 to 1.1 2.22 × 10−4 Eq. (6) 275

m3 s−1 MPa W ◦ C−1

LSHE

Tair,in Tair,out UA

22 36 40

◦

Superheat

5

◦

IHE

C C W ◦ C−1

◦

C

Fig. 5. Time evolution of the building load.

A. Nguyen, P. Eslami-Nejad and M. Badache et al. / Energy & Buildings 202 (2019) 109343

7

Since every point of the compression cycle is calculated at each time step, one can evaluate the energy consumption of the compressor and the auxiliary system at any given time. The transient performance of the DX-GSHP system is therefore evaluated in terms of the system’s overall COP (heat pumps and auxiliary system), calculated with:



t

COP (t ) = 0 0

(QHP (t ) + QAUX (t ) )dt t

(27)

(W (t ) + QAUX (t ) )dt

where QHP is the load provided by the heat pump, QAUX is the load provided by electrical auxiliary system and W is the work of the compressor.

Fig. 7. Time evolution of the ground temperature at z = 1.2 m.

4. Results and discussion Fig. 6 (a) illustrates the evolution of the DX-GSHP system overall COP (Eq. (27)) while Fig. 6(b) presents the cumulative energy consumption of the system (i.e. compressor and auxiliary heater together), during the 4-month heating period for all four simulated cases. The energy consumption profiles are calculated using the denominator in Eq. (27). The slope of these profiles varies continuously due to the change of the heating demand (Fig. 5) and the operating conditions of the DX-GSHP. The total system and auxiliary system energy consumption calculated at the end of the 4-month simulation for all the mentioned cases are tabulated in Table 3. Note that the auxiliary heating of all the simulated cases are comparable, and accounts for a small portion of the load provided. Of all cases, the first two (with IHE) present the lowest energy con-

Fig. 8. Vertical profile at March 1st for the a) fluid b) and borehole wall temperature. The arrows indicate the direction of the flow.

Fig. 6. a) Time evolution of the average system COP and b) the total system energy consumption. Table 3 Transient performance of the simulated cases.

Case Case Case Case

#1 #2 #3 #4

Total consumption (kWh)

Auxiliary consumption (kWh)

Overall COP (-)

1586.4 1576.4 1927.2 1787.3

86.8 84.2.8 118.2 93.1

3.66 3.68 3.01 3.25

sumptions and highest COPs, thus highlighting the positive effect of the IHE on the efficiency of a DX-GSHP. As for the effect GSTV, it is always detrimental to the performance of the DX-GSHP, with different extent however. When comparing case #1 and #2, both performance curves are similar throughout the entire simulation. The overall performance is only slightly lower for case #1 (with GSTV). However, when comparing cases #3 and #4, performance diverges after the 1200th h. This is due to the cold front penetration from the ground surface as the winter progresses (Fig. 7). The calculated average COP of Case #3, which accounts for GSTV, is 7% lower than Case #4. This means the lower shallow ground temperature reduces the heat extraction capacity of the BHE. For better understanding, detailed analysis of the temperature profile in the BHE is presented and described as follows. Fig. 8 (a) and (b) present respectively the fluid and borehole wall temperature profiles for March 1st. For cases case #3 and #4 (without an IHE), fluid temperature increases respectively from −15.5 ◦ C and −12.5 ◦ C to reach the ground temperature (10 ◦ C), then, decreases in the second leg to exit the borehole at respectively 11 ◦ C and −8.5 ◦ C. This shows that thermal short-circuiting with the first leg as well is significant. This very curvy shape is due to the fluid’s transition from two-phase to single-phase. As most heat is extracted from the ground during the evaporation (first 15 m of the first leg), the heat exchange between the ground

8

A. Nguyen, P. Eslami-Nejad and M. Badache et al. / Energy & Buildings 202 (2019) 109343

Fig. 9. Vertical profile at March 1st for the vapor quality. The arrows indicate the direction of the flow.

and the single-phase fluid is relatively marginal and is merely to provide superheated vapor to the compressor. The cold soil at the upper part of the borehole in case #3 makes this even more difficult. Therefore, this operation (particularly at low flow) makes the BHE to operate inefficiently. On the other hand, for cases #1 and #2 (with an IHE), the fluid temperature profile is quasi-constant, meaning that evaporation is occurring over the entire length of the U-pipe and heat extraction from top to bottom of the borehole is almost uniform. Note that the small temperature increase is attributed with the geothermal gradient as well as pressure gain due to gravity. Interestingly, GSTV has little effect on fluid and borehole wall temperature profiles with an IHE is employed. Among all the cases from #1 to #4, the lowest fluid temperature is obtained in cases #3, as the system needs to provide superheated vapor at the oulet of the borehole in a colder environment because of GSTV.

Fig. 9 shows that the vapor quality profiles are different with and without an IHE. These profiles provide a better understanding of the evaporation process occurring in the BHE. It indicates that, for cases #3 and #4 (without an IHE), CO2 enters the borehole at a relatively high inlet vapour quality (x = 0.45), fully evaporates at 12 m, superheats to a quality of 1.1 at the bottom of the BHE and finally, loses heat to the first leg to exit with a 5 ◦ C superheat (x = 1.02). For cases #1 and #2, the IHE allows heat to be recovered from the low pressure side, thereby allowing the fluid to enter the borehole at a much lower quality (x = 0.03) and subsequently to steadily increase to x = 0.7. This confirms a quasi-uniform heat extraction rate over the borehole length in the presence of an IHE. Again, one can see that GSTV has little impact on the operation of the heat pump when IHE is used. Fig. 10 presents a visual appreciation of the ground temperature contours (isotherms) on March 1st for all four cases. Naturally, more heat extraction from the ground means lower ground temperature and farther cold front penetration from the borehole wall. For instance, Fig. 10(d) shows that, although the seasonal ground temperature variation is not accounted for, the ground temperature is colder in the first 20 m. However, in Fig. 10(b) the isotherms are subvertical, which is due to uniform heat extraction from top to bottom of the borehole when the IHE is employed. A similar trend is observed for cases where GSTV is taken into account (cases #1 and #3). Fig. 11 presents the heat pump compression cycle using the CO2 pressure-enthalpy diagram on March 1st for all four cases. The yellow curve represents the 22 ◦ C isotherm while the orange curve indicates the 5 ◦ C superheat. At a glance, all cases release heat at almost similar discharge pressure below the critical point at 7–7.2 MPa (subcritical). This is due to the following reasons; 1) heat is delivered at a relatively low temperature (36 ◦ C) and 2) the

Fig. 10. Close up of the temperature profile at March 1st for: a) Case 1, b) Case 2, c) Case 3 and d) Case 4.

A. Nguyen, P. Eslami-Nejad and M. Badache et al. / Energy & Buildings 202 (2019) 109343

9

approach temperature and, thereby increase the system COP. Furthermore, the small discrepancy with Eq. (6) is due to numerical tolerance for faster convergence. 5. Conclusion

Fig. 11. Compression cycle on the pressure-enthalpy diagram for CO2 at March 1st. The yellow curve represents the 22 ◦ C isotherm while the orange curve indicates the 5 ◦ C superheat. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 12. Distribution of P2 for Case 1.

building load on March 1st is relatively low (less than one third of the peak). As shown in Fig. 11, case #1 and case #2 have identical compression cycles, which confirms the negligible effect of GSTV on system operation when the heat pump is equipped with the IHE. This is proven by calculating the cycle COP; 4.74 and 4.76 for case #1 and case #2 respectively. The heat pump cycle presents significantly lower evaporating temperatures in cases #3 and #4 (without an IHE) compared to first two cases. As mentioned above, it is due to the drop of CO2 pressure at the borehole inlet to provide the required superheat at the outlet. Among cases #3 and #4, the difference between heat pump cycles is also notable. This confirms that GSTV negatively affects heat pump operations when an IHE is not used. The calculated cycle COP is 2.89 and 3.12 for case #3 and case #4, respectively. These results clearly show that equipping the heat pump with an IHE effectively counters the effect of GSTV and thermal short-circuiting in the BHE. In other words, adding an IHE means that the BHE no longer needs to provide the required superheat, thus, improving its efficiency. Another benefit of using an IHE is a pressure drop reduction in the BHE, as displayed in Fig. 11. In terms of the discharge pressure control, Fig. 12 presents a scatter plot (case #1) of P2 with respect to T3 for the whole simulation period. This shows that the DX-GSHP operated in both subcritical and transcritical cycle during the simulation. One should note that T3 ranged from 27.5 to 31 ◦ C, while the return air temperature was fixed at 22 ◦ C. These high approach temperatures may be attributed to the refrigerant mass flow, discharge pressure and pinch point location in the LSHE, as suggested by [6]. Increasing the inlet refrigerant pressure (and temperature) could decrease the

The ground temperature is susceptible to affect the performance of a DX-GSHP due to the required amount of superheat at the suction of the compressor. In this study, the effect of an IHE on the performance of a CO2 DX-GSHP operating in heating mode was assessed through transient energy simulations. To do so, the DXGSHP model developed by [20] was extended to account for an IHE and GSTV. First, the numerical approach for modeling GSTV was validated by comparing the proposed model against a reference solution. Results showed very good agreements (MAE = 0.02 ◦ C) for the computation of the ground temperature profile under the effect of GSTV. Then, a 4-month heating transient simulations were performed for four cases: with and without and IHE, and GSTV. It was shown that equipping the heat pump with an IHE could mitigate thermal short-circuiting by maintaining a two-phase flow in both legs of the BHE, thereby keeping a relatively constant heat flux in the BHE. This was true for cases with and without the effect of GSTV. In terms of performance, the overall system COP was 22% higher in the case with an IHE when GSTV was considered (3.66 versus 3.01), a significant gain for a simple addition. This finding is contradictory to the one in [10], showing thereby the importance of more complex models and transient simulations. As for the modeling, it was shown that, for short-term simulations (4month here), the effect of GSTV is marginal when an IHE is used. Otherwise, providing superheated vapor to compressor was shown to be more difficult, as the ground temperature close to the surface is relatively colder. For long-term simulations, care should be taken before neglecting the effect of GSTV, even if the heat pump is equipped with an IHE, as the BHE is relatively short. Moreover, discharge pressure was controlled using the correlation proposed by [25]. For this study, results indicate that increasing the discharge pressure could further increase the system COP by reducing the approach temperature in the LSHE. These finding open the door for future research on controls for DX-GSHP in dynamic transient analysis. The developments presented in this work represent a step in this direction. Naturally, experimental validations are recommended to support the conclusions presented here. Finally, the proposed model can be used to simulate the transient operation of a CO2 DX-GSHP and evaluate the performance of systems under various controls, ground conditions and building loads. Declaration of Competing Interest We wish to confirm that there are no known conflicts of interest associated with this publication and there as been no significant financial support for this work that could have influenced its outcome. Acknowledgments The authors would like to acknowledge the funding received by the Office of Energy Research and Development (OERD) (grant no. EU-BE-22) of Canada through the Energy Innovation Program (EIP) to undertake this research. References [1] B.T. Austin, K. Sumathy, Parametric study on the performance of a direct-expansion geothermal heat pump using carbon dioxide, Appl. Therm. Eng. 31 (2011) 3774–3782. [2] M. Badache, P. Eslami-Nejad, M. Ouzzane, Z. Aidoun, L. Lamarche, A new modeling approach for improved ground temperature profile determination, Renew. Energy 85 (2016) 436–444.

10

A. Nguyen, P. Eslami-Nejad and M. Badache et al. / Energy & Buildings 202 (2019) 109343

[3] D. Bauer, W. Heidemann, H.J. Diersch, Transient 3D analysis of borehole heat exchanger modeling, Geothermics 40 (4) (2011) 250–260. [4] I.H. Bell, J. Wronski, S. Quoilin, V. Lemort, Pure and pseudo-pure fluid thermophysical property evaluation and the open-source thermophysical property library coolprop, Ind. Eng. Chem. Res. 53 (2014) 2498–2508. [5] A. Bidarmaghz, G.A. Narsilio, I.W. Johnston, S. Colls, The importance of surface air temperature fluctuations on long-term performance of vertical ground heat exchangers, Geomech. Energy Environ. 6 (2016) 35–44. [6] Y.G. Chen, Pinch point analysis and design considerations of CO2 gas cooler for heat pump water heaters, Int. J. Refrig. 69 (2016) 136–146. [7] M. Cimmino, P. Eslami-Nejad, A simulation model for solar assisted shallow ground heat exchangers in series arrangement, Energy Build. 157 (2017) 227–246. [8] E.A. Elias, R. Cichota, H.H. Torriani, Q. de Jong van Lier, Analytical soil-temperature model, Soil Sci. Soc. Am. J. 68 (2004) 784–788. [9] P. Eslami-Nejad, M. Ouzzane, Z. Aidoun, A quasi-transient model of a transcritical carbon dioxide direct-expansion ground source heat pump for space and water heating, Appl. Therm. Eng. 91 (2015) 259–269. [10] H. Ghazizade-Ahsaee, M. Ameri, Effects of using expander and internal heat exchanger on carbon dioxide direct-expansion geothermal heat pump, Appl. Therm. Eng. 136 (2018) 389–407. [11] H. Ghazizade-Ahsaee, M. Ameri, Energy and exergy investigation of a carbon dioxide direct-expansion geothermal heat pump, Appl. Therm. Eng. 129 (2018) 165–178. [12] Y. Guo, G. Zhang, J. Zhou, J. Wu, W. Shen, A techno-economic comparison of a direct expansion ground-source and a secondary loop ground-coupled heat pump system for cooling in a residential building, Appl. Therm. Eng. 35 (2012) 29–39. [13] L. Jensen-Page, G.A. Narsilio, A. Bidarmaghz, I.W. Johnston, Investigation of the effect of seasonal variation in ground temperature on thermal response tests, Renew. Energy 125 (2018) 609–619. [14] T.H. Lim, R.D. De Kleine, G.A. Keoleian, Energy use and carbon reduction potentials from residential ground source heat pumps considering spatial and economic barriers, Energy Build. 128 (2016) 287–304. [15] V.J. Lunardini, Heat Transfer in Cold Climates, Van Nostrand Reinhold Co, Toronto, Canada, 1981. [16] D. Marcotte, P. Pasquier, F. Sheriff, M. Bernier, The importance of axial effects for borehole design of geothermal heat-pump systems, Renew. Energy 35 (2010) 763–770.

[17] G. Mihalakakou, On estimating soil surface temperature profiles, Energy Build. 34 (2002) 251–259. [18] D. Ndiaye, Reliability and performance of direct-expansion ground-coupled heat pump systems: issues and possible solutions, Renew. Sustain. Energy Rev. 66 (2016) 802–814. [19] A. Nguyen, P. Eslami-Nejad, M. Badache, A. Bastani, Pressure-enthalpy coupled thermal resistance and capacity model (PH-TRCM) for direct-expansion borehole heat exchangers: application for supercritical CO2 , Geothermics 76 (2018) 50–59. [20] A. Nguyen, P. Eslami-Nejad, A transient coupled model of a variable speed transcritical CO2 direct expansion ground source heat pump for space heating and cooling, Renew. Energy 140 (2019) 1012–1021. [21] A. Nguyen, P. Pasquier, D. Marcotte, Influence of groundwater flow in fractured aquifers on standing column wells performance, Geothermics 58 (2015) 39–48. [22] A. Nguyen, P. Pasquier, D. Marcotte, Thermal resistance and capacity model for standing column wells operating under a bleed control, Renew. Energy 76 (2015) 743–756. [23] A. Nguyen, P. Pasquier, D. Marcotte, Borehole thermal energy storage systems under the influence of groundwater flow and time-varying surface temperature, Geothermics 66 (2017). 110–8 [24] P. Pasquier, D. Marcotte, Joint use of quasi-3D response model and spectral method to simulate borehole heat exchanger, Geothermics 51 (2014) 281– 299. [25] P.C. Qi, Y.L. He, X.L. Wang, X.Z. Meng, Experimental investigation of the optimal heat rejection pressure for a transcritical CO2 heat pump water heater, Appl. Therm. Eng. 56 (2013) 120–125. [26] I. Sarbu, C. Sebarchievici, General review of ground-source heat pump systems for heating and cooling of buildings, Energy Build. 70 (2014) 441–454. [27] F. Tang, H. Nowamooz, Long-term performance of a shallow borehole heat exchanger installed in a geothermal field of alsace region, Renew. Energy 128 (2018) 210–222. [28] (2017). The MathWorks inc. Matlab version 9.3.0 (R2017b). Natick. Massachusetts. [29] W. Wu, H.M. Skye, Progress in ground-source heat pumps using natural refrigerants, Int. J. Refrig. 92 (2018) 70–85. [30] A. Zarrella, M. Scarpa, M. De Carli, Short time step analysis of vertical ground– coupled heat exchangers: the approach of caRM, Renew. Energy 36 (9) (2011) 2357–2367.