Transient simulation and multi-objective optimization of a VSD ground source heat pump in various usage

Energy Conversion and Management 197 (2019) 111847

Contents lists available at ScienceDirect

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

Transient simulation and multi-objective optimization of a VSD ground source heat pump in various usage

T

⁎

Soheil Kaviana, Cyrus Aghanajafib, , Nader Dizadjia a b

Department of Mechanical Engineering, Science and Research Branch, Islamic Azad University, Tehran, Iran Faculty of Mechanical Engineering, K. N. Toosi University of Technology, Tehran, Iran

A R T I C LE I N FO

A B S T R A C T

Keywords: MOPSO Exergy destruction Irreversibility Variable speed drive Ground source heat pump LCC

The main purpose of this research is to optimize the life cycle cost (LCC) while minimizing the annual irreversibility of the system in three different modes of a vertical ground heat pump (water-to-water) including heating, cooling, and multi-usage. The multi-objective optimization of the LCC and exergy destruction was done by using the multi-objective particle swarm optimization (MOPSO) algorithms and also by considering all involving effective parameters such as the number of boreholes, depth of buried pipes in ground source heat exchanger (GSHX), fan coil setpoint temperature and surface area of condenser and evaporator. In order to reduce energy consumption and meet the building's heat requirement at any time, the compressor was simulated with variable speed drive (VSD) in transient operations. In this study, the investigated refrigerant was R134a, and the shell-and-tube type of heat exchanger was considered for the condenser and evaporator. The obtained results show that the use of the heat pump in single heating mode leads to lower exergy destruction and better performance of the system in terms of exergy efficiency. Also, the economic results indicate that the total initial costs for the multi-usage mode is 27% higher than the single heating mode while it is 8.5% higher than the single cooling mode. The average electricity cost per operation time in multi-usage is 15% higher than single heating and 3% lower than single cooling and equal to 0.139 $/hour.

1. Introduction One of the largest consumers of energy is known to be the residential buildings sector, where heating and cooling comprise a major portion of their energy consumption [1]. In recent years, the use of clean energy for residental buildings has been gaining considerable attention due to the problems stemming from fossil fuel pollutants. Furthermore, rising fuel prices highlight the importance of renewable energy sources as clean and sustainable alternatives to fossil fuels [2]. The country's investment in geothermal energy almost doubled from 2010 to 2015, and the use of geothermal energy has increased by 45% [3]. One of the common ways to extract this clean energy source for heating and cooling of residential buildings is ground source heat pumps (GSHP) [4]. GSHP systems use the ground as a heat source/sink medium to meet the demand for heating or cooling load [5]. Energy efficiency and initial investment cost are among the most critical points in the design of geothermal heat pumps. Hence, thermoeconomic optimization and study of various parameters of heat pump design play essential roles in this regard. Accordingly, Hepbasli et al. [6] presented an energy and exergy

⁎

analysis of a GSHP system in heating mode with a 50-m vertical U-bend ground heat exchanger and the case study was evaluated for a 65-m2 room in the Solar Energy Institute, Ege University, Izmir in Turkey. In this regard, relations for the energy and exergy analysis of a closed loop vertical GSHP system for different components were derived using average measured parameters obtained from experimental results and mass, energy, exergy, and entropy balance equations. Sayyaadi Et al. [7] presented a model based on steady-state simulation and thermodynamic and thermoeconomic optimization using an evolutionary algorithm for a vertical ground-source heat pump system in cooling mode with a constant refrigerant flow rate. Their approach was applied to minimize the total economic cost based on the total revenue requirement (TRR) method and the exergy destruction of the system. Three levels of optimization, namely single thermodynamic objective, single economic objective, and multi-objective optimizations, were applied and the results showed that the levelized fuel cost for the multi-objective, thermoeconomic, and thermodynamic optimized systems are respectively 50%, 22.5% and 54.4% less than the base case. Liu et al. [8] simulated an air-source variable refrigerant flow (VRF) system and GSHP that used single-stage scroll compressors and vertical ground loop

Corresponding author. E-mail addresses: [email protected], [email protected] (S. Kavian), aghanajafi@kntu.ac.ir (C. Aghanajafi), [email protected] (N. Dizadji).

https://doi.org/10.1016/j.enconman.2019.111847 Received 2 May 2019; Received in revised form 16 July 2019; Accepted 17 July 2019 0196-8904/ © 2019 Elsevier Ltd. All rights reserved.

Energy Conversion and Management 197 (2019) 111847

S. Kavian, et al.

Nomenclature A B Bo C C1 C2 Cinitial Cp Ctube CL CTP d DS e E Ė

f fo & m fsalv G h H i İ k L LCC ṁ Nt Nu NTU P Pij Pr PR PWF Q̇ r R Re s Ṡ T t vij VG W Ẇ x Xij Z

Greek symbols

Surface area (m2) Baffle spacing (m) Boiling number Heat capacity (W/K) Personal learning coefficient Social learning factor Initial cost of the system ($) Specific heat capacity (kJ/(kg K)) the clearance between adjacent tubes (m) Tube layout constant Tube count calculation constant Diameter (m) Shell internal diameter (m) Discount rate Two-phase convection multiplier The rate of exergy (kW/h) Friction factor Operating and maintenance fraction Salvage fraction Mass flux (kg/(m2s)) Heat transfer coefficient (W/(m2 K)) Enthalpy (kJ/(kg K)) Interest rate Rate of irreversibility (kW) Thermal conductivity (W/(mK)) Length (m) Life cycle cost ($) Mass flow rate (kg/s) Number of tubes Nusselt number Number of transfer units Pressure (kPa) Best position Prandtl number Pitch ratio Present Worth Factor Heat transfer rate (W) uniform random numbers Thermal resistance (m2K/W) Reynolds number Specific entropy (kJ/kg K) Rate of entropy (kW/K) Temperature (K) Thickness (m) Velocity of the swarms Compressor inlet volume (m3) Inertia weight Power (kW) Vapour quality Position of the swarms Investment cost ($)

ηelec ηis ηmech ηtotal ηR ηv μ ε ρ ω

Compressor electrical efficiency Isentropic efficiency Compressor mechanical efficiency Compressor total efficiency Exergy efficiency Volumetric efficiency Viscosity (kg/(ms)) Heat exchanger effectiveness Density (kg/m3) Compressor speed (rpm)

Subscripts 0 1 2 3 4 5 6 7 8 annual b C comp cond dest eff eq evap ExV f fanc g gen gshx h HP in isen K l Max min p out ref sheat TP v w

Environmental state (dead state) Evaporator outlet Condenser inlet Condenser outlet Evaporator inlet Ground heat exchanger inlet Ground heat exchanger outlet Fan coil inlet Fan coil outlet Annual energy consumption bore Cold fluid flow Compressor Condenser Destroyed Heating system efficiency Equivalent Evaporator Expansion valve Fouling resistance Fan coil best global position Generation Ground heat exchanger Hot fluid flow Heat pump Inlet Isentropic Location Liquid Maximum Minimum pipe Outlet Refrigerant Space heating Two-phase Vapour Water

operational costs associated with exergy destruction. Their results showed that the electricity consumption of the ITES system was 9.8% lower in comparison with a conventional system. Hakkaki-Fard et al. [5] evaluated the energy efficiency, system costs and the relative payback period of an air-source heat pump (ASHP) compared with directexpansion ground-source heat pump (DX-GSHP) in a residential building in the cold climate city of Montreal. Their results showed that the energy consumption of the DX-GSHP system could be reduced by 50% and the relative payback period of the GSHP is more than 15 years

heat exchanger in heating and cooling modes. They also compared the systems in terms of energy efficiency for a small office building between two different US climates and reported that the typical GSHP has better performance and, in general, GSHP is more energy-efficient than the air-source VRF system, especially when the building has significant heating loads. On thermo-economic optimization, Sanaye et al. [9] evaluated a thermo-economic optimization of an ice thermal energy storage system (ITES) for air-conditioning applications using genetic algorithm, in which the objective function included the capital and 2

Energy Conversion and Management 197 (2019) 111847

S. Kavian, et al.

longer than ASHP. Hoseini Rahdar et al. [10] studied a model based on optimization of R-717 and R-134a ice thermal energy storage air conditioning systems using non-dominated sorting genetic algorithm II (NSGA-II) and multi-objective particle swarm optimization (MOPSO) to obtain the optimal design parameters leading to the optimal objective functions, exergy efficiency and total cost rate. Their results showed that by using multi-objective optimization of parameters, the electricity consumption of the R-717 refrigerant system decreases by 11%, with 8% more reduction than the R-134a refrigerant system. Liu et al. [11] evaluated the feasibility and performance of the hybrid ground-source heat pump system (HGSHP) for an office building in a cold area and three systems, including one HGSHP system with boiler as auxiliary heat source, ground source heat pump system (GSHP), and a traditional electric water chiller. The transient simulation of the system in hourly dynamic heating and cooling modes was modeled in the TRNSYS software. Their results indicated that the HGSHP system with the boiler as auxiliary heat source operates more efficiently and reduces traditional energy consumption. Zhang et al. [12] performed exergy destruction and thermo-economic assessment of a novel combined heat and power (CHP) system integrating partial biomass gasification with GSHP system in heating mode, and the Aspen Plus software was used in their simulation. Their results showed that the exergy efficiency of the system reached about 13.65% and the unit cost of hot water, the power of gas turbine and the steam turbine were obtained to be 92.5 $/GJ, 12.4 $/GJ, and 23.9 $/GJ, respectively. Li et al. [13] presented a model based on multi-objective optimization of the previously proposed integration system by genetic algorithm, with consideration of energy, economic and environmental performances. There were three optimal parameters: the rated power of the power generation unit, carbon conversion ratio and the set temperature of hot water. The results showed that the primary energy saving ratio, annual total cost saving ratio, CO2 emission reduction ratio and performance indicator are 7.61%, 23.62%, 66.52%, and 32.58%, respectively. Chahartaghi et al. [14] studied energy and exergy and thermo-economic optimization in a GSHP with closed-loop horizontal heat exchanger for domestic water heating. Modeling was performed in steady state by applying engineering equation solver (EES) software and the effect of variations in inlet water temperature of the evaporator for several types of refrigerants on the coefficient of performance (COP) of the system was investigated. Their results showed that R507a had the highest COP in comparison with the other refrigerants. Also, R507a had the minimum total annual cost (TAC) and the maximum cost of the system belonged to the ground heat exchanger, which equaled 34% of the total cost. Patel et al. [15] evaluated the thermodynamic and ecological performance of a reversed-Brayton-cycle-operated heat pump in a steady state through many-objective and multi-objective optimization by the heat transfer search (HTS) algorithm. Maximization of the COP, ecological coefficient of performance, exergy efficiency and ecological function of the heat pump were considered as objectives and four operating variables of the heat pump, namely temperature ratio, compressor pressure ratio, and hot side and cold side heat exchanger effectiveness, were investigated in their optimization study. Different decision-making approaches were used to select a final optimal solution from the Pareto optimal set of many-objective optimization, and finally the results revealed that the deviation between optimization and experimental results of the heat pump system are 10.95% and 12.3% in COP and exergy efficiency, respectively. Despite all these developments and studies on the GSHP, there is no comprehensive study on the multi-objective optimization of the irreversibility and economic performance of the variable-speed drive heat pump systems in transient simulation and comparison of heating, cooling, and double-use applications. In this paper, the transient simulation of a variable speed drive (VSD) heat pump in three different modes—heating, cooling, and multiusage— is the main focus. Here, the evaluation of the three systems is prioritized over mere comparison and the primary purpose of analyzing

Fig. 1. The schematic of the primary heat pump cycle with ground source energy and fan coil unit equipment in heating mode.

heating, cooling, and multi-usage modes is to understand how the values of essential variables could change in these types of systems, so the designers could decide what the best system is for their needs from the point of view of climate region and economic investment. The main purpose of this research is to optimize the Life Cycle Cost (LCC) while minimizing the annual irreversibility of the system. The general multiobjective optimization of the LCC and exergy destruction was done according to the MOPSO algorithms. The optimization was conducted considering all involving effective parameters including the number of boreholes, depth of buried pipes in GSHX, fan coil setpoint temperature, condenser and evaporator surface area. 2. Methodology The purpose of this study is the evaluation of a heating and cooling system with ground source energy and water-to-water heat pump for a residential building, as shown in Fig. 1. A U-Tube model heat exchanger was used for the ground source, while for both condenser and evaporator, shell-and-tube configuration was utilized. The refrigerant used in condenser and evaporator was selected to be R-134a and the secondary fluid in the GSHX and fan coil units was water. The compressor was of the scroll-type with a variable speed drive (VSD) able to synchronize with desired load for a fan coil unit at each time step. In order to simulate the heat pump behavior, a numerical code was developed in the Engineering Equation Solver (EES) software [16] and for evaluation of transient response of ground source heat exchanger and the total system, TRNSYS software was used. For this purpose, a new type in TRNSYS software [17] was developed by C++ programming that could communicate with input and output data of the EES program at different times during the simulation. The general multiobjective optimization of the LCC and exergy destruction was done according to the MOPSO algorithms and with the aid of the MATLAB program. For this purpose, a code was written in MATLAB that ran TRNSYS software in each simulation. 2.1. Case study and hourly loads simulation The case study in this research is a residential building in the city of Tehran, Iran with three floors and seven thermal zones. The total area of the thermal zones is 192 m2, and the building specifications are presented in Table 1. The weather information for this region was obtained from the climate meteorological data. Hourly and average ambient temperature changes are presented in Fig. 2. The minimum and maximum ambient temperatures are 4 °C in January and 42 °C in July, respectively. Average monthly relative humidity (RH) (%) and total radiation on the horizontal ground surface (kWh/(m2.month)) are presented in Fig. 3. 3

Energy Conversion and Management 197 (2019) 111847

S. Kavian, et al.

Table 1 Building Specifications. Material

Orientation

Thickness/ Property

U-value (W/ m2.k)

Total Area (m2)

Exterior Walls

NorthEast NorthWest SouthEast SouthWest Horizontal Horizontal NorthEast NorthWest SouthEast SouthWest

0.305 m 0.305 m 0.305 m 0.305 m 0.445 m 0.7 m Double layer Double layer Double layer Double layer

1.291 1.291 1.291 1.291 0.973 0.704 2.83 2.83 2.83 2.83

83.2 98.52 61.97 75.42 170 55.32 35.4 6 1.8 34.1

Roof Floor Windows

Fig. 3. Monthly average RH and total radiation of Tehran.

The hourly heating and cooling loads of the building were simulated with these data as input in the TRNBuild software [18].

n

Qtotal =

(1)

i=0

2.2. Heat pump modeling

Outlet water and refrigerant temperature are obtained from the NTU method by the following relations [19–21]:

The heating cycle consists of a closed loop system in which the water absorbs heat from the refrigerant in the condenser and then transfers the heat to the air in the fan coil. In the cooling cycle, the evaporator absorbs the heat of the space through fan coil units. A numerical model was developed for simulating heat pump performance with R-134a refrigerant, and the following assumptions were made:

• Saturated liquid and saturated vapor conditions are assumed at the • • • • • •

∑ Qi

̇ Q̇ = εQmax

(2)

Q̇max = Cmin (Th, in − Tc, in )

(3)

̇ p C = mc

(4)

where Th and Tc are the hot and cold sides of the inlet fluid, respectively, ε represents the effectiveness of the heat exchanger and can be calculated by the following equation in a single phase refrigerant and shell-and-tube type heat exchanger [20,22]:

condenser outlet and compressor inlet (point 3 and 1, repsectively in Fig. 1) respectively. The expansion valve undergoes an isenthalpic process (point 3 to 4). The compressor undergoes an isentropic process (point 1 to 2). Pressure drop in connecting tubes is neglected. Flow for the refrigerant inside tubes and water inside the shell is one-dimensional (for both condenser and evaporator). Heat loss to the surroundings is assumed to be negligible. The directions of heat transfer to the system and work transfer from the system are positive.

−1

1 + exp[−N (1 + C 2)1/2] ⎫ ε = 2 ⎧1 + C + (1 + C 2)1/2 × ⎨ 1 − exp[−N (1 + C 2)1/2] ⎬ ⎩ ⎭

(5)

For a two-phase condensation or boiling process, it can be assumed that the fluid temperature stays essentially constant, or in other words, the fluid acts as if it had infinite specific heat and all the heat-exchanger effectiveness relations are expressed as [19–22]:

ε = 1 − e−NTU 2.2.1. Shell-and-tube heat exchanger Shell-and-tube heat exchanger was assumed for both the condenser and the evaporator. A section-by-section one-dimensional scheme for calculating local heat transfer in the heat exchanger was used. The schematic layout of the control volume of the heat exchanger is presented in Fig. 4. Total heat transfer of the heat exchanger is the sum of heat transfer of each section, which can be calculated by the following equations:

N = NTU =

UA Cmin

(6)

(7)

In the mentioned formula, U and A are the overall heat transfer coefficient and heat exchanger surface area, respectively and can be calculated as [23]:

Aout = πdout LNt

Fig. 2. Hourly and monthly average ambient temperature of Tehran. 4

Energy Conversion and Management 197 (2019) 111847

S. Kavian, et al.

Fig. 4. The schematic of the control volume of the heat exchanger.

Q̇ = εQ̇max Q̇ = εQ̇max

(8)

obtained from the correlations listed in Table 3.

1 1 = Rt + U hout

(9)

2.2.2. Compressor A variable speed drive with a scroll compressor is considered in this study. The refrigerant mass flux can be obtained by the following equations [21,29,30]:

1 d t d + Rf , in ⎤ out + wall out Rt = Rf , out + ⎡ ⎥ ⎢ h d k Dm in in wall ⎦ ⎣ Dm =

(10)

ṁ ref = ηv × ρ1 × VG × ω

dout − din ln(dout / din )

(18)

(11) 2

P P ηv = 0.9207 − 0.0756 ⎛ 2 ⎞ + 0.0018 ⎛ 2 ⎞ ⎝ P1 ⎠ ⎝ P1 ⎠

where do and di are outside and inside tube diameters, Nt is the number of the tubes, ho and hi are shell-and-tube heat transfer coefficients, tw and kw are thickness and conductivity of the wall, and Rfi and Rfo are inside and outside fouling resistance values. The main characteristics of the condenser and evaporator heat exchanger are listed in Table 2. The specifications of the shell-and-tube heat exchangers and the length to diameter ratio (L/Ds) are chosen according to the suggestion by Jiang et al. [24] and commercial manufacturing. For the shell side, the McAdams correlation for heat transfer coefficient is used which can be obtained as [23]:

DG hout = 0.36 ⎜⎛ e s ⎟⎞ ⎝ μ ⎠

0.55

1/3

⎛ Cp μ ⎞ ⎝ k ⎠

⎜

⎟

ṁ Gs = As As =

DS CB PT

Ctube = PT − dout Ds2 CTP ⎞ Nt = 0.785 ⎛ 2 ⎝ CL ⎠ (PR)2dout

⎛ μb ⎟⎞ ⎝ μw ⎠

⎜

0.14

⎛k⎞ ⎝ De ⎠

⎜

⎜

2 4(PT2 − πdout /4) πdout

⎜

⎟

(19)

in which, ω is not constant and varies with time and heating loads. As the compressor is assumed to undergo an isentropic process, the power can be calculated as [7,21]:

̇ Wcomp =

ṁ ref (h2, isen − h1 ) ηtotal

(20)

ηtotal = ηis ηmech ηelec

(21)

(12)

P ηis = 0.85 − 0.046667 ⎛ 2 ⎞ ⎝ P1 ⎠

(13)

So the coefficient of performance (COP) of the heat pump can be obtained by the following relations:

(14)

COPHP =

⎟

⎜

⎟

(22)

̇ Qcond ̇ Wcomp

(23)

(15) Table 2 Main Characteristics of shell and tube heat exchangers.

(16)

where Ds is the internal shell diameter, PT is the pitch size, and B is the baffle cut (baffle spacing), which is considered 25%, CTP is the tube count calculation constant which is considered 0.90 for two tube passes, CL is the tube layout constant which is considered 1.0 for the 90° layout [23], and PR is pitch ratio (PT/dout). De is the equivalent diameter and for the square pitch can be calculated as [23]:

Deq =

⎟

(17)

The single and two-phase Nusselt numbers inside the tube are 5

Item

condenser

evaporator

Number of tube passes Number of shell passes Inner tube diameter (m) Outside tube diameter (m) Tube layout (o) Baffle cut (%) Inside fouling resistance (m2.k/W) Outside fouling resistance (m2.k/W) Tube wall conductivity (W/m.k) The ratio of length to diameter

2 1 0.0082 0.0095 90 25 0.00018 0.00009 111 8

2 1 0.0082 0.0095 90 25 0.00018 0.00009 111 8

Energy Conversion and Management 197 (2019) 111847

S. Kavian, et al.

the final step, the compressor drive speed (ω ) was calculated using Eq. (18) and having the refrigerant mass flow rate as a known parameter.

Table 3 Heat transfer Correlations inside the tube. Author

Type and phase

Correlation

Gnielinski [21,23]

The single phase inside the tube

Nu =

The Boiling flow inside the tube

f= h = Ehl hl = 0.023Rel0.8 Prl0.4 kl/ din

Gungor and Winterton [25–27]

COPSystem =

One of the approaches to the economics process is to use the life cycle cost method that takes into account all future expenses. This method provides a means of comparison of future costs with today’s costs. As a result, the present energy cost in periods can be calculated by the following relation [32,33]:

Rel = G (1 − x ) din/ μl

E = 1 + 3000Bo0.86 + 1.12 Cavallini and Zecchin [23,28]

2.3. Life cycle cost method

(f / 2)(Re − 1000) Pr 2 1 1.07 + 12.7(f / 2) 2 ⎜⎛Pr 3 − 1⎟⎞ ⎝ ⎠ (1.58 ln Re − 3.28)−2

The Condensation flow inside the tube

0.75 ρ x ⎛⎜ l ⎞⎟ 1−x ρ ⎝ g⎠

( )

0.41

PWF (N , i, e ) =

0.8 0.33 hTP = 0.05Reeq Prl kl/ din μ

ρ

Reeq = Rev ⎛ v ⎞ ⎛ l ⎞ ⎝ μl ⎠ ⎝ ρv ⎠

0.5

+ Rel

̇ Qcond ̇ Wcomp + Ẇ pumps + Ẇ fans

1 ⎡ 1 + i ⎞N ⎤ 1−⎛ (e − i) ⎢ ⎝1 + e ⎠ ⎥ ⎦ ⎣

(25)

where e is the market discount rate, i is the energy inflation rate, and N is periods in years. Finally, the life-cycle cost (LCC) can be obtained as [34]:

LCC

(24)

1 + i ⎞N ⎤ ⎡ = Cinitial ⎡1 + fo & m × PWF − fsalv ⎛ +⎢ ⎢ ⎝1 + e ⎠ ⎥ ⎣ ⎦ ⎣

2.2.3. Ground source heat exchanger For evaluation of transient and local behavior of the ground and heat exchanger with time, the borehole was simulated with type 557 TRNSYS software [17]. The details of the borehole are shown in Table 4. The numbers of ground layers were considered to be four, with their properties shown in Table 5 [31].

P ⎛ ̇ ⎞⎤ ⎜ (Qannual ) × η ×PWF⎟ ⎥ eff ⎝ ⎠⎦

(26)

where Cinitial is the initial cost of the system, fo & m is the operating and maintenance fraction, fsalv is the salvage fraction, p is the energy price, ̇ is the annual energy ηeff is the heating system efficiency, and Qannual consumption. In the above equation, the initial costs of the components are obtained from the equations [9,35,36] in Table 6.

2.2.4. Model simulation The heat pump simulation was performed using a code developed in the EES software [16]. The required input and output data are as below. The required input data of the program are:

2.4. Exergy analysis

• Inlet water temperature to the condenser (T ) • Inlet water temperature to the evaporator (T ) • Inlet water mass flow rate to the condenser • Inlet water mass flow rate to the evaporator • The desired heating load • Number of tubes in condenser and evaporator • Length of the condenser and evaporator • Refrigerant Type 8

In a closed-loop system, irreversibility or exergy destroyed is an important parameter because, the lower the irreversibility, the higher the efficiency. The general exergy destroyed or the irreversibility rate can be calculated as [37–39]:

6

̇ I ̇ = EẊ , dest = T0 Sgen

(27)

The rate of entropy generation can be obtained by the following relation [37]:

The calculated output data of the program are:

• • • • • •

̇ = Sgen

Outlet water temperature from the condenser (T7) Outlet water temperature from the evaporator (T5) The refrigerant mass flow rate and drive speed COP of the heat pump Compressor power exergy destroyed of the heat pump

∑ ṁ out sout − ∑ ṁ in sin − ∑

Qk̇ Tk

(28)

where Qk̇ is the thermal load entering the system and Tk is the temperature of its location. Hence, the irreversibility rate of each component of the system can be obtained from the following relations [6] in Table 7. The exergy efficiency of the whole system can be defined as the ratio of the desired exergy output to the exergy input [6,7,40], which can be determined as:

For evaluation of transient behavior of the GSHX and the total system, TRNSYS software was used. For this purpose, a C++ code was developed to facilitate communication between the output data of the EES program and the input data of TRNSYS and vice versa in different runs. For this purpose, with the assist of C++ coding, a new type in TRNSYS software is created that writes input data into a text file, then runs the EES software and finally reads the output data from another text file. In EES software, a code was developed that reads input data from the text file created by TRNSYS software and writes the output data into another text file that can be read by the typed created in TRNSYS. The calculation steps and the general algorithm are presented in Fig. 5. The refrigerant mass flow rate and other output data were obtained by having known input parameters like desired heating load and calculating in iterations loop cycles, as shown in the flow chart. In

Table 4 The specifications of the borehole heat exchanger.

6

Item

Value

Borehole radius (m) The outer radius of the U-tube pipe (m) The inner radius of the U-tube pipe (m) Reference Borehole flow rate (kg/hr) Pipe thermal conductivity (W/m.k) Fill thermal conductivity (W/m.k) Center-to-center half distance (m) Fluid specific heat (kJ/kg.K)

0.0625 0.0194 0.016 612.36 0.43 1.6 0.032 4.19

Energy Conversion and Management 197 (2019) 111847

S. Kavian, et al.

Table 5 The properties of the ground in different layers. Layer

Depth (m)

Thermal conductivity (W/m.K)

Volumetric heat capacity (kJ/m3.K)

Clay Limestone (upper carbonate layer) Limestone (lower carbonate layer) Limestone (upper carbonate layer)

0–16 16–30 30–61 61-∞

1.832 3.11 2.6 1.92

2550 4198 4198 4198

ηR =

EẊ , desired EẊ , used

Table 6 Component costs function.

(29)

Component

T0 ⎞ ̇ EẊ , desired = ⎜⎛1 − ⎟ Qsheat Tin, air ⎠ ⎝

Pump

(30)

T ̇ EẊ , used = Wcomp + Ẇ fan + Ẇ pump + ⎛1 − 0 ⎞ Q̇gshx T soil ⎠ ⎝ ⎜

Price function ($) 0.71 Zpump = 705.48 × Ẇ pump ⎛1 + ⎝ 0.6123 Zevap = 16648.3 × Aevap ⎜

Evaporator

⎟

(31)

2.5. Multi-objective particle swarm optimization

Expansion valve

ZExV = 114.5 × ṁ ref

condenser Compressor

Zcond = (516.621 × Acond ) + 268.45

GSHX

As a multi-objective extension of PSO, MOPSO gives a set of Pareto solutions [41]. In this method, an arrangement of solutions is considered and characterized by their position and velocity vectors. The best solution for each particle is stored in memory as an experience, and the best solution among all of the particles is achieved and called the best global particle [10,41]. In any iteration of the algorithm, the position and velocity of the swarms are updated based on the following

0.2 ⎞⎟ 1 − ηpump

Zcomp =

39.5ṁ ref 0.9 − ηis

⎠

( ) ln ( ) P2 P1

P2 P1

Zgshx = 1.1Lp + 19.1Lb Lp = 2Lb

equations [10,41]:

vij = Wvij + C1 r1 (Pij − Xij ) + C2 r2 (gi − Xij )

(32)

Xij = Xij + vij

(33)

Fig. 5. Logic flow chart of the heat pump cycle. 7

Energy Conversion and Management 197 (2019) 111847

S. Kavian, et al.

Table 7 The exergy destroyed of the component. Component

The irreversibility or the exergy destroyed

Compressor

̇ Icomp = T0 ṁ ref (s2 − s1)

Condenser

̇ Icond = T0 [ṁ ref (s3 − s2) + ṁ w (s7 − s8)]

Evaporator

̇ Ievap = T0 [ṁ ref (s1 − s4 ) + ṁ w (s5 − s6)]

Expansion valve

̇ IExV = T0 ṁ ref (s4 − s3)

Fan coil units

̇ I fanc = T0 ⎡ṁ w (s8 − s7) + ⎣

GSHX

̇ Igshx = T0 ⎡ṁ w (s6 − s5) − ⎣

Q̇fanc ⎤ Tin, air ⎦ Q̇gshx Tsoil

⎤ ⎦

Table 8 MOPSO parameters for modeling.

Fig. 7. Comparison of the boiling heat transfer coefficient predictions with experimental measurements [27].

Parameter

Value

Number of iterations Population size Inertia weight Global learning coefficient Personal learning coefficient

50 15 0.7298 1.491 1.491

Table 9 Range of changes for selected variables. Variable

Minimum Value

Maximum Value

Number of boreholes Borehole depth (m) Fan coil heating setpoint (C) Fan coil cooling setpoint (C) Tube length (Condenser) (m) Tube length (Evaporator) (m)

1 20 35 5 1 1

10 300 80 19 10 10

Fig. 8. Comparison of the condensation heat transfer coefficient of current study with the predicted results of the Cavallini et al. [45] model.

where, W is the inertia weight, C1 is the personal learning coefficient, C2 is the social learning factor, r1 and r2 are the uniform random numbers, Pij is the best position and gi is the best global position. The setting parameters of MOPSO for this modeling are specified, at an appropriate range, in Table 8 [10,42]. Minimizations of LCC and annual exergy destruction are considered as the objective functions for the optimization. Six parameters, namely the number of boreholes, depth of buried pipes of the GSHX, fan coil setpoint temperature in cooling and heating mode, and tube lengths of condenser and evaporator are considered as decision variables, and the range of changes are presented in Table 9. The other unknown geometric parameters of the heat exchanger like the number of tubes and surface area of the heat exchanger can be obtained by the aid of Eq. (8),

(15) and (16) and having known parameters from optimization (tube length) and the aspect ratio of heat exchanger (L/Ds) from Table 2. There have been some studies that used the coupling of TRNSYS with other programs that contain multi-objective optimization algorithms such as TRNSYS with MATLAB and GenOpt [42,43], or TRNSYS with GenOpt [44]. In this study, a code in MATLAB was developed according to the algorithm presented in Fig. 6. For this purpose, with the assist of C++ programming, a new TRNSYS type was created that writes the objective functions into a text file. In MATLAB, several subprograms were developed for reading text output of TRNSYS data, printing the value of the selected optimization variables into the .dck TRNSYS file, MOPSO algorithm, and running the TRNSYS software.

Fig. 6. Schematic flowchart of the coupling of the MATLAB and TRNSYS programs for the MOPSO method. 8

Energy Conversion and Management 197 (2019) 111847

S. Kavian, et al.

Table 10 The properties and experimental condition of [21]. Component

Volumetric flow (m3/h)

Total number of tubes

Number of tube passes

Tube length (m)

Inner tube diameter (m)

Outside tube diameter (m)

condenser evaporator

1.5 1.14

20 76

2 2

0.8 0.8182

0.013 0.00822

0.016 0.00952

Fig. 9. Comparison of COP of the current study with experimental measurements [21].

Fig. 11. Pareto front: best optimum values for the objective functions (the total irreversibility and LCC).

3. Model validation

et al. [21]. The experimental measurements were based on a water-towater heat pump with shell-and-tube heat exchangers and variable speed compressor drive. In order to compare the current model with experimental results, the input data were set the same as the experimental condition as is shown in Table 10. Also in this simulation, the inlet temperature of the evaporator and compressor speed were considered 302 K and 575 rpm, respectively. The COP results at different condenser inlet temperatures are presented in Fig. 9. The average accuracy for the different condenser temperatures between the predicted COP of the simulated data and the experimental measurement is about 7.72%, which is acceptable assuming neglibible heat loss in the pipe.

The boiling heat transfer coefficient predictions for R134a refrigerant were compared with experimental measurements of Saitoh et al. [27] inside a 1.12-mm-diameter tube for the mass flux of 150 kg/ (m2s) and the heat flux of 15 kW/m2 (Fig. 7). As it is shown, the predicted results agree well with an average deviation of 6.69% compared with the experimental measurements. The condensing heat transfer coefficient of R134a refrigerant was compared with the predicted results of the Cavallini et al. [45] model inside an 8-mm-diameter tube for the mass flux of 200 kg/(m2s) at 40 °C saturation temperature. Also, the difference between saturation and wall temperature was considered constant and equal to 4 °C. Their model predicted values with an average deviation of 5% with respect to experimental results. The results of the comparison of the current study predictions with those of Cavallini et al. [45] are presented in Fig. 8. As shown, the predicted heat transfer coefficients in this study agree well with an average deviation of 5.66%. In order to verify the theoretical model of the heat pump, the simulation results of the EES code for R134a refrigerant were compared with the data from the experimental study of the Mendoza-Miranda

4. Results and discussion The results are divided into the heating and cooling load demands, multi-objective optimization results, thermodynamic analysis, comparison of the performance of different operation modes and economic analysis. In Section 4.1, hourly simulation results of the heating and cooling load of the building are presented. In Section 4.2, the results of the multi-objective optimization with the purpose of minimizing the

Fig. 10. Hourly simulation of the heating and cooling loads of the building. 9

Energy Conversion and Management 197 (2019) 111847

S. Kavian, et al.

Table 11 The results of optimized values for different operation modes. Variables

Multi-usage (heating & cooling)

Heating mode

Cooling mode

Number of boreholes Borehole depth (m) Fan coil heating setpoint (C) Fan coil cooling setpoint (C) Tube length (Condenser) (m) Number of tubes (Condenser) Tube length (Evaporator) (m) Number of tubes (Evaporator) Average Exergetic efficiency Total annual exergy destructions (kW) LCC ($)

7 50.96 47.53 9.79 1.57 192 1.47 169 0.105 17150.95 35846.18

5 60.00 47.07 – 1.40 153 1.08 90 0.168 6473.00 20400.82

6 51.76 – 12.26 1.36 143 1.57 192 0.064 8746.73 27895.90

17.66 kW, taking place in July. Different operation modes of the heat pump were simulated by the aid of type 2b of TRNSYS software as a control signal for switching the heat pump on or off. In heating mode, comfort setpoint temperature of the building was considered to be 21 °C; hence, below this temperature, the heat pump would be turned on. In cooling mode, comfort setpoint temperature of the building was considered 23 °C, above which the heat pump would be turned on.

4.2. Multi-objective optimization Because of transient simulation of this model for one year, one of the challenges of this study was a long computation time for multiobjective optimization even with a few iterations; the simulation was thus performed with a high-performance server with 2 × Intel® Xeon® Gold 6130 16-Core CPU and 64 GB DDR4 ECC RAM for faster calculations. The Pareto optimum frontier chart in multi-objective optimization for the three different modes of heat pump operation (heating, cooling and multi-usage) is presented in Fig. 11. Selection of the best optimum points exist on Pareto front needs; a decision-making process depends on the importance of each objective for designers. In other words, it is possible to define a weighting coefficient for each objective according to its importance. In this study, it was assumed that both objectives had the same value and the weighting coefficient was considered to be one for each of them. Selection of the best point was done by the suggested method of Sayyaadi et al. [7] with the aid of introducing a hypothetical yellow point in Fig. 11. This is an ideal point at which both objectives have their optimum value and therefore, this point does not exist practically. Therefore, the closest point of the Pareto frontier to the yellow point that has the shortest distance might be considered as a

Fig. 12. Exergy destructions of the components for different operation modes.

irreversibility and LCC for various usage of the heat pump are discussed. In Section 4.3, exergy destruction for different components and performance of different conditions are investigated. Finally, in Section 4.4, initial costs, operation costs, and LCC in 20 years of the cooling, heating, and multi-usage of the heat pump are discussed. 4.1. Heating and cooling demand Hourly simulation of heating and cooling loads of total residential building is shown in Fig. 10. The results were calculated using hourly temperature, solar radiation, and RH. In this simulation, the signs of heating and cooling loads were considered positive and negative, respectively. According to the diagram, the heating peak load is 13.22 kW and occurs in January. Also, the peak load of the cooling mode is

Fig. 13. The monthly average value of COP and compressor speed in different months for the modes under study. 10

Energy Conversion and Management 197 (2019) 111847

S. Kavian, et al.

Fig. 14. Hourly heating peak load simulation of the performance of the heat pump at 13th of January.

Fig. 15. Hourly cooling peak load simulation of the performance of the heat pump at 22th of July.

desirable optimum solution which is shown in green in the figure. The results of decision variables for optimum points are presented in Table 11. The obtained results show that the total annual irreversibility values of the multi-usage mode are 62% and 49% higher than the corresponding values in heating and cooling modes respectively. In the optimized operation mode, the LCC values of 20 years of multi-usage are 43% and 22% higher than the corresponding values in heating and cooling modes, respectively.

Table 12 The result of the initial costs of different components for various cases. Component

Multi-usage

Heating mode

Cooling mode

Compressor price ($) Condenser price ($) Evaporator price ($) Expansion valve price ($) GSHX price ($) Fan coil price ($) Fan coil pump price ($) GSHX pump price ($)

268.3875 2464.108 1267.237 19.53364 7598.583 1831.433 180.7208 146.7749

271.8433 1793.971 564.7896 19.67358 6390 1019.262 145.2831 127.4088

300.3533 1634.994 1505.358 20.19133 6615.184 2453.918 222.4442 132.7707

4.3. Thermodynamics analysis The mean annual exergy destruction of each component in terms of operation times of the heat pump is presented in Fig. 12. The results 11

Energy Conversion and Management 197 (2019) 111847

S. Kavian, et al.

4.4. Economic analysis The life-cycle period of simulation was considered 20 years and calculations were done based on an inflation rate of 17% and an interest rate of 13% according to the Iranian national data. Electricity cost for the power of components including the compressor, fan, and pumps was considered to be 0.06 $/kWh with an inflation rate of 7% for each year. The results of the initial costs of different components for various cases are presented in Table 12. According to the results, it can be concluded that the highest initial cost belongs to the GSHX component and also the prices decrease by 15% and 12% for single usage in heating and cooling modes in comparison to the multi-usage mode. Fig. 16 shows the initial and annual electricity costs for the different operation modes. Annual energy costs for multi-usage, at 1211$, is 66% higher than single heating mode and 41% higher than single cooling mode. Also, the total initial costs for multi-usage is 27% higher than the single heating mode and 8.5% higher than the single cooling mode. Since the heat pump in the multi-usage mode operates over the entire year, it is evident that more energy is consumed. Consequently, the LCC and annual electricity costs will be higher than the heating and cooling mode systems, but the benefits for the building residents are much more in the multi-usage mode because they will profit from thermal comfort during the entire year. For a fair comparison of electricity costs of three types of system usages, the annual electricity costs were divided into operation times, and the results are shown in Table 13. According to the results, the average electricity costs per operation time ($/hour) in multi-usage is 15% higher than single heating mode and 3% lower than single cooling mode.

Fig. 16. Comparison of initial and annual electricity costs for different cases. Table 13 The average electricity costs per operation time of the heat pump for different usages. Operation

Electricity costs ($/hour)

Multi-usage Cooling Heating

0.139 0.143 0.117

show that the highest irreversibility belongs to the condenser component in heating mode. In cooling and multi-usage modes, the fan coil unit has the highest exergy destruction. Exergy destruction of the condenser in the heating mode is 82% higher than the cooling mode. The irreversibility of the GSHX in cooling mode is 70% higher than heating mode and 19% higher than multi-usage. It can also be concluded that the GSHX is more efficient in heating compared with cooling. For the evaporator, exergy destruction in the cooling mode is 87% higher than the heating mode. Also for the compressor, irreversibility in the heating mode is 32% higher than cooling, which is due to the high inlet temperature of the fan coil units in the heating loads. Fig. 13 presents the monthly average amounts of the COP and compressor speed of the simulated heat pump for different cases including heating, cooling, and multi-usage. The results show that the average compressor speed is higher in January and July for heating and cooling modes, respectively. In January, the optimum average COP values for the heat pump are 4.49 and 4.31 for multi-usage and heating modes, respectively. The average speeds of the compressor are 5346 and 5596 rpm for multi-usage and heating modes, respectively. In July, the optimum average COP values for the heat pump are 4.46 and 4.10 for multi-usage and cooling modes, respectively. Also, the average speeds of the compressor are 5176 and 5276 rpm for multi-usage and heating modes, respectively. The average annual COP values for multiusage, heating and cooling modes are 5.58, 4.67 and 4.71 respectively. Fig. 14 shows the hourly performance of the heat pump at the heating peak load on 13 January. The peak load at 8 AM occurs with the amount of 12.9 kW. At this time, COP values decrease to 3.75 and 3.57 for multi-usage and heating modes, respectively, as expected, while the pressure ratios increase to 4.65 and 4.89, respectively. Also, the outlet temperatures of the GSHX or the inlet temperature of evaporator decrease to 12.1 °C and 11.1 °C for multi-usage and heating modes, respectively. Fig. 15 shows the hourly performance of the heat pump at the cooling peak load on 22 July. The peak load At 3 PM occurs with the amount of 17.58 kW. At this time, the COP decreases to 2.83 and 2.34 for multi-usage and cooling mode as expected and also the pressure ratio increases to 4.75 and 5.18 respectively. Also, the outlet temperature of the GSHX or the inlet temperature of condenser increases to 35.1 °C and 38.5 °C for multi-usage and cooling modes, respectively.

5. Conclusions In this research, a multi-objective optimized approach for achieving the best design was adopted by investigating single and multi-usage of the vertical ground source heat pump in cooling and heating modes. The optimization process was carried out using a MOPSO algorithm. The multi-objective optimization focused on limited energy and monetary resources, simultaneously. The proposed method covered both thermodynamic and economic aspects of the system design and the component selection. Irreversibility and LCC for the ground source heat pump components were determined based on derived equations in a transient state and also for meeting the demand load at each time step, the driver speed of the compressor was considered to be variable. The main conclusions are summarized in the following points:

• For the heating, the optimum setpoint temperature of the fan coil was about 47 °C. • For the cooling, based on usage type (multi or single usage) the • • • • • • 12

optimum setpoint temperature of the fan coil was between 9 °C and 12 °C. The optimized number of boreholes was between 5 and 7 with a depth of 50–60 m. Because of more needed boreholes in multi-usage mode, the costs of GSHX increased by up to 15% and 12% compared with single heating and cooling modes, respectevely. The obtained optimum results of average exergy efficiency for multiusage, heating, and cooling were 10.5%,16.8%, and 6.4% respectively. From the point of view of exergy efficiency, the use of the heat pump in single heating mode results in lower exergy destruction and better performance of the system. According to the optimized results, the average annual COP were 5.58, 4.67 and 4.71 for multi-usage, heating and cooling modes, respectively. Since the heat pump in the multi-usage mode operates over the entire year, it is evident that more energy is consumed, and consequently in multi-usage mode, the 20-year LCC values were 43% and

Energy Conversion and Management 197 (2019) 111847

S. Kavian, et al.

22% higher than the heating and cooling modes, respectively.

[17] Klein S, et al. TRNSYS 16–A TRaNsient system simulation program, user manual. Solar Energy Laboratory. Madison: University of Wisconsin-Madison; 2004. [18] SEL T, CSTB T, Multizone Building modeling with Type56 and TRNBuild; 2004, TRNSYS. [19] Browne M, Bansal P. An elemental NTU-ε model for vapour-compression liquid chillers. Int J Refrig 2001;24(7):612–27. [20] Holman J. Heat transfer. 10th ed. Mc-GrawHill Higher education; 2010. [21] Mendoza-Miranda JM, et al. Variable speed liquid chiller drop-in modeling for predicting energy performance of R1234yf as low-GWP refrigerant. Int J Refrig 2018;93:144–58. [22] Mills A. Heat Transfer. Richard Drwin Inc; 1992. [23] Kakac S, Liu H, Pramuanjaroenkij A. Heat exchangers: selection, rating, and thermal design. CRC Press; 2002. [24] Jiang F, Deng X. Deep heat transfer performance and ratio of length to diameter in shell and tube heat exchangers. International Conference on Information and Management Engineering. Springer; 2011. [25] Gungor K. Simplified general correlation for saturrated flow boiling and comparison of correlation with data. Chem Eng Res Dev 1987;65:148. [26] Gungor KE, Winterton R. A general correlation for flow boiling in tubes and annuli. Int J Heat Mass Transf 1986;29(3):351–8. [27] Saitoh S, Daiguji H, Hihara E. Correlation for boiling heat transfer of R-134a in horizontal tubes including effect of tube diameter. Int J Heat Mass Transf 2007;50(25–26):5215–25. [28] Browne M, Bansal P. An overview of condensation heat transfer on horizontal tube bundles. Appl Therm Eng 1999;19(6):565–94. [29] Brown JS, Yana-Motta SF, Domanski PA. Comparitive analysis of an automotive air conditioning systems operating with CO2 and R134a. Int J Refrig 2002;25(1):19–32. [30] Oritz T, Li D, Groll EA. Evaluation of the performance potential of CO2 as a refrigerant in air-to-air Air conditioners and heat pumps: system modeling and analysis. ARTI; 2003. Final Report. [31] Olfman MZ, Woodbury AD, Bartley J. Effects of depth and material property variations on the ground temperature response to heating by a deep vertical ground heat exchanger in purely conductive media. Geothermics 2014;51:9–30. [32] Duffie JA, Beckman WA. Solar engineering of thermal processes. John Wiley & Sons; 2013. [33] Kavian S, Saffari Pour M, Hakkaki-Fard A. Optimized design of the district heating system by considering the techno-economic aspects and future weather projection. Energies 2019;12(9):1733. [34] Jaber S, Ajib S. Optimum, technical and energy efficiency design of residential building in Mediterranean region. Energy Build 2011;43(8):1829–34. [35] Kavanaugh SP. Simulation and experimental verification of vertical ground-coupled heat pump systems. Oklahoma: State University; 1985. [36] Valero A, et al. CGAM problem: definition and conventional solution. Energy 1994;19(3):279–86. [37] Sonntag RE, et al. Fundamentals of thermodynamics Vol. 6. New York: Wiley; 1998. [38] Shakib SE, et al. Multi-objective optimization of a cogeneration plant for supplying given amount of power and fresh water. Desalination 2012;286:225–34. [39] Razmi A, et al. Thermodynamic and economic investigation of a novel integration of the absorption-recompression refrigeration system with compressed air energy storage (CAES). Energy Convers Manage 2019;187:262–73. [40] Habibi M, Hakkaki-Fard A. Long-term energy and exergy analysis of heat pumps with different types of ground and air heat exchangers. Int J Refrig 2019;100:414–33. [41] Coello CC, Lechuga MS. MOPSO: a proposal for multiple objective particle swarm optimization. Proceedings of the 2002 Congress on Evolutionary Computation. CEC'02 (Cat. No. 02TH8600). IEEE; 2002. [42] Rey A, Zmeureanu R. Multi-objective optimization of a residential solar thermal combisystem. Sol Energy 2016;139:622–32. [43] Asadi E, et al. A multi-objective optimization model for building retrofit strategies using TRNSYS simulations, GenOpt and MATLAB. Build Environ 2012;56:370–8. [44] Magnier L, Haghighat F. Multiobjective optimization of building design using TRNSYS simulations, genetic algorithm, and Artificial Neural Network. Build Environ 2010;45(3):739–46. [45] Cavallini A, et al. Condensation in horizontal smooth tubes: a new heat transfer model for heat exchanger design. Heat Transfer Eng 2006;27(8):31–8.

• Annual electricity costs for multi-usage was 66% higher than single heating mode and 41% higher than the single cooling mode. • The average electricity costs per operation time ($/hour) in multi-

usage was 15% higher than single heating mode and 3% lower than single cooling mode.

Therefore, it can be concluded that the proper design and optimization can greatly increase the system efficiency for a residential building, but in general, the application of the ground source heat pump can be more beneficial in cold climates for meeting the heating loads from in terms of economic costs and irreversibility of the system. Declaration of Competing Interest None. References [1] Mosleh HJ, Hakkaki-Fard A, DaqiqShirazi M. A year-round dynamic simulation of a solar combined, ejector cooling, heating and power generation system. Appl Therm Eng 2019;153:1–14. [2] Mosleh HJ, et al. A new desalination system using a combination of heat pipe, evacuated tube and parabolic trough collector. Energy Convers Manage 2015;99:141–50. [3] Lund JW, Boyd TL. Direct utilization of geothermal energy 2015 worldwide review. Geothermics 2016;60:66–93. [4] Habibi M, Hakkaki-Fard A. Evaluation and improvement of the thermal performance of different types of horizontal ground heat exchangers based on technoeconomic analysis. Energy Convers Manage 2018;171:1177–92. [5] Hakkaki-Fard A, et al. A techno-economic comparison of a direct expansion groundsource and an air-source heat pump system in Canadian cold climates. Energy 2015;87:49–59. [6] Hepbasli A, Akdemir O. Energy and exergy analysis of a ground source (geothermal) heat pump system. Energy Convers Manage 2004;45(5):737–53. [7] Sayyaadi H, Amlashi EH, Amidpour M. Multi-objective optimization of a vertical ground source heat pump using evolutionary algorithm. Energy Convers Manage 2009;50(8):2035–46. [8] Liu X, Hong T. Comparison of energy efficiency between variable refrigerant flow systems and ground source heat pump systems. Energy Build 2010;42(5):584–9. [9] Sanaye S, Shirazi A. Thermo-economic optimization of an ice thermal energy storage system for air-conditioning applications. Energy Build 2013;60:100–9. [10] Rahdar MH, et al. Modeling and optimization of R-717 and R-134a ice thermal energy storage air conditioning systems using NSGA-II and MOPSO algorithms. Appl Therm Eng 2016;96:217–27. [11] Liu Z, et al. Feasibility and performance study of the hybrid ground-source heat pump system for one office building in Chinese heating dominated areas. Renewable Energy 2017;101:1131–40. [12] Zhang X, et al. Exergetic and exergoeconomic assessment of a novel CHP system integrating biomass partial gasification with ground source heat pump. Energy Convers Manage 2018;156:666–79. [13] Zhang X, et al. Optimization analysis of a novel combined heating and power system based on biomass partial gasification and ground source heat pump. Energy Convers Manage 2018;163:355–70. [14] Chahartaghi M, et al. Energy and exergy analyses and thermo-economic optimization of geothermal heat pump for domestic water heating. Int J Low-Carbon Technol 2019. [15] Patel VK, Raja BD. A comparative performance evaluation of the reversed Brayton cycle operated heat pump based on thermo-ecological criteria through many and multi objective approaches. Energy Convers Manage 2019;183:252–65. [16] Klein SA, Alvarado F. EES: Engineering equation solver for the Microsoft Windows operating system. 1992. F-Chart software.

13