Sizing and performance analyses of a combined heating and cooling system with the integration of short- and long-term storages

Technological advancements in the utilization of renewable energy sources have unveiled potentials for increasing building energy efficiency. Integrating heat pump-based energy systems with thermal storages is a suitable option to meet the thermal requirements of modern buildings and exploiting the available renewable energy sources. However, how to size the main components of a heat pump-based energy system with the integration of short- and long-term storages is not yet well explored. Therefore, this study focused on the design and performance analyses of an integrated heating and cooling system consist of a heat pump, borehole long-term thermal storage, and hot water tank short-term thermal. Heat pump models were introduced as parametric models based on the producer data. The dynamic thermal model of the energy system was developed and analysed in MATLAB. Different combinations of heating and cooling loads were tested. Integration of cooling and heating systems was discussed through different operation strategies and challenges were addressed. The results of the parametric analysis identified the key parameters affecting the design of components and efficiency of the system. Moreover, the results showed that lower cooling to heating load ratio leads to an excessive reduction of the ground temperature and overall efficiency over the long-term operation.


Introduction
Energy consumption in building sector accounts for approximately 40% of total energy use and one third of the direct and indirect CO 2 emissions in the world [1,2] [3]. Cooling demand in many modern buildings is now increasing [4]. The heating and cooling energy may be simultaneously required. Hence, the combination of thermal loads is another aspect to be considered when designing the building energy systems [5,6].
Heat pumps-based technologies can be coupled with available renewable energy sources such as geothermal energy to simultaneously cover heating and cooling demands [7,8]. Due to stable ground temperature below certain depths, borehole thermal energy storage (BTES) can store the heating energy in the ground for long periods despite the low specific heat of the ground material [9,10]. BTES often yields slow charging or discharging rates. Coupling BTES with short term storage units could smoothen the short-term operations [5,9,11].
There is a considerable amount of research on design and performance analysis of heating and cooling energy systems [7,[12][13][14][15][16][17][18][19]. There are well-established methods for sizing the elements of a complex energy system stand alone. However, how to size the main components of a heat pump-based energy system with the integration of short-and long-term storages is not yet well explored. In this regard, a general methodology for the preliminary design of heat pump-based building energy systems is highly demanded. In this study, an integrated heating and cooling system is described and analysed. This analysis aimed to investigate the most influencing parameters in sizing and overall efficiency improvement of the energy system. This paper is structured as follows. The description of the energy system and its operation modes is given in Section 2. In Section 3 methods for modelling of the system parts are explained. Results are presented in Section 4, followed by a discussion in Section 5, and the main conclusions of the research in Section 6.

System description
The system considered in this study consisted of a heat pump, substations, a hot water storage tank (HWST), and a BTES. The heating and cooling loads were assembled into one lumped demand. Heat sources for the heat pump could be the heat from the cooling system or BTES. The heating energy supplied by the heat pump was primarily stored in HWST. Surplus heat was stored in the BTES via a heat exchanger. The building heating system received the required heating energy from the HWST. The minimum supply temperature for the heating system was set to 50°C. A district heating connection was also considered as an auxiliary heat supply unit for the peak load coverage and if required for necessary temperature lift. A schematic of the system is shown in Figure 1. Hourly values for annual heating and cooling loads were calculated or measured in advance and they were used as the model inputs. The model was tested with two different sets of cooling and heating loads with distinct characteristics, named Load 1 and Load 2 to study the transferability of the model outcomes. Data for the Load 1 was calibrated based on the thermal demands of barracks of a military base, while data for Load 2 was collected from the simulation of a hospital energy system in Norway. Hourly values of thermal loads are presented in Figure 2 and Figure 3.  Load 1 had relatively high and constant heating and cooling demand throughout the year, as shown in Figure  2. Load 2 had lower total annual energy demand and the heating and cooling demands were not overlapping, as shown in Figure 3. Table 1 gives key information for the users. The ratio of the total cooling demand to heating demand, α, was introduced to include the combination of loads into the assessment. For the base case, Load 1 was characterized by higher peak heating load and a high α value, while Load 2 had a larger peak cooling load and lower α value. The energy system had three operating modes depending on the heating and cooling demands. Heating mode: The system operated in heating mode when the heating demand was higher than the cooling demand. The heat pump extracted energy from the ground. The BTES was used as the heat source for the heat pump. The actual evaporator load was calculated based on the borehole extraction heat rate. The return brine from the evaporator was firstly sent to the cooling system and then to the BTES. The condenser heat was sent to the HWST.
Cooling mode: The system operated in the cooling mode when the cooling demand was higher than the heating demand. The heat pump was connected to the cooling system as the heat source. Therefore, in the case, the actual evaporator load was considered based on cooling demand. The condenser heat rate in this mode was higher than the heating demand. Depending on the energy level of the HWST, there could be a surplus heat on the condenser side of the system. In that case, the ground acted as a heat sink, where the excess condenser heat was delivered to the BTES.
Free cooling: The system operated in free cooling mode when the heating demand was zero. The heat pump in this mode was off. The only flow of energy was from the cooling system to the BTES.

Methods
In this section, the approaches chosen for the sizing of each component of the system is described in detail. The dynamic thermal balance model of the system developed in MATLAB is explained. Further, a sensitivity analysis of the key parameters of the system is presented.

Heat pump model
Heat pump models may be obtained from parametric modelling of the refrigerant cycle. On the other hand, empirical heat pump models are developed from the records of the capacities and compressor power in relation to one or several operating conditions [20][21][22]. The heat pump model used for this study was based on the regression of simulation data generated from the compressor manufacturer BITZER. Nominal capacities were derived as functions of evaporation temperature , compressor discharge temperature , and nominal cooling capacity , as follows: For each of the target users, a heat pump was selected with the specifications as listed in Table 2. The evaporation temperature and the compressor discharge temperature were defined as: and ∆ were the assumed minimum temperature differences at the secondary side of the condenser and the evaporator, respectively. The inlet temperatures on the secondary sides were calculated in the system model. The outlet temperatures on the secondary sides were set as an input to the heat pump model.
The part load operation of the heat pump was calculated by considering actual loads on evaporator and condenser. Actual evaporator load was calculated according to the operation mode of the system. Actual condenser load was determined based on the real-time storing capacity of the HWST. The mass flow rate of the water circulating the condenser side was determined as: where ∆ ℎ is the estimated 40K temperature drop between the condenser inlet and the outlet of the condenser. Cooling and heating part loads were defined as ratios of the actual evaporator and the condenser loads to their corresponding nominal capacities. To consider the simultaneity between the cooling and the heating part loads, the highest part load was chosen as the overall part load via equations 8-10.
There were no limitations set regarding step control or minimum speed control for the part load, hence it was assumed that the compressor could operate at any part load between zero and one. The overall coefficient of performance (COP) for heating and cooling modes was calculated as: Beside the COP of the heat pump stand-alone, the total coefficient of performance of the system , was also examined to represent the overall annual efficiency of the heat pump-based energy system. Hence, the annual energy input from the peak load and the storage units were considered as:

Borehole thermal energy storage model
Research on heat transport inside and outside the boreholes have revealed the importance of thermal response of ground heat exchangers in design procedures [10,17,19,23,24]. In addition, adequate knowledge of local geology and thermal properties of the ground is demanded to optimize the storage process under various conditions [1,3,5,8,9]. For sizing the borehole configuration, an Energy Earth Designer (EED) evaluation was carried out to assess the thermal properties of the ground and heat carrier fluid in relation to the thermal loads and sizing dimensions [16]. The ground properties were assumed from the typical values in Norway [25,26]. The thermal properties of the ground are presented in Table 3. The sizing of the BTES was performed by considering the nominal capacity of evaporator and the heat extraction rate per length from the ground to calculate the total length of boreholes as: The number of the boreholes was then decided by dividing the total length of effective boreholes with the chosen length for the boreholes. The BTES configuration was modelled as a control volume using a deterministic dynamic approach. The main parts of the BTES were single U-tube pipes, the borehole with filling material, and the surrounding ground. To simplify the geometry of the borehole field, a square configuration was considered. The thermal energy balance used to model the BTES was given as:

Hot water storage tank model
The role of HWST in the present study was to enable the supply of heating energy independent of heat pump operation. Sizing approaches of the hot water storage tanks addressed in literature are diverse. Depending on the usage profile and composition of energy system different methods can be adopted [11,27]. The volume of the HWST was calculated based on the maximum energy coverage of the HWST as: where ℎ , is the design heat rate. Design heat rate was chosen to correspond to the largest heat demand for 90% of the hours of a year. The design heat rate for Load 1 and Load 2 was 60% and 90% of maximum heat demand, respectively.
The HWST was modelled by dividing the water in the tanks into two identical horizontal sections. These sections were assumed to have a uniform temperature with the inlet and outlet of the tanks at the tank top and bottom. Internal heat exchange between the sections was neglected. Thermal stratification was simplified by assuming a maximum 40K temperature difference between two sections and linearly reducing when the temperature at the upper section falls below 60°C [28,29]. The thermal energy balance in the tank was calculated by summation of the inlet and the outlet flows and heat losses to ambient as:

Substation models
The substation model for the heating system in this study was developed by considering heat supply from hot water tank and district heating acting as peak load covering unit as: The peak load kicked in when the temperature of the water from the hot water tank drops below 50°C or when the tank heat rate is not covering the required heat rate. A maximum of 40K temperature drop in the heating system was assumed to estimate the mass flow rate in the heating system. Substation model for the cooling system was similarly developed by considering maximum 5K temperature gain from the cooling system to calculate the flow rate in the cooling system.
The models of the components were connected to establish the dynamic model of integrated heating and cooling system. Heating and cooling loads, as well as the heat pump operating temperature range and performance coefficients, were provided as input for the simulations. Distribution heat losses and hydraulic components were not considered in the model.

Sizing and performance analysis
The response of the sizing and performance outputs of the system model to changes in the design parameters was studied by One-At-a-Time sensitivity analysis [30][31][32]. Sizing variables considered to investigate were volume, total effective length, and maximum heat extraction rate of the BTES as well as energy coverage rate of the HWST. The performance indicators considered were and , . For the sensitivity analysis, 20% increase and decrease were applied to the parameters listed in Table 4. The changes were applied to one parameter at a time. The only exception was the heat pump performance coefficients that were changed for 10%, because of the valid capacity range of the chosen heat pumps.

Results
Result of this study is presented in three parts. First, the model described in previous sections was tested on both Load 1 and Load 2 to reveal the main annual performance characteristics of the system. Then the result of sensitivity analysis is presented. Finally, the performance analysis of the 10-years operation of the system is presented.

Annual performance analysis
The initial temperature of the BTES and the HWST was set to 8.5°C and 40°C, respectively. Figure 4 shows the hourly thermal performance of the BTES in one-year operation. In general, the temperature oscillations of Load 1 were less than Load 2. This could be due to the presence of cooling demand throughout the year. In the heating mode hours, BTES temperature was decreasing due to heat extraction from the ground. In the cooling mode operation with Load 1, the BTES was charged with surplus heat from the condenser. Due to generally high heat demands, as well as the large HWST volume and undersized heat pump, BTES charging in the cooling mode operation was not enough to recover the temperature to the initial level. Cooling demands of Load 2 were supplied only in free cooling mode. Hence, the charging rate and consequently temperature uplift of the BTES was higher for Load 2. Heating mode operation of Load 2 was characterized as the heating only condition. Therefore, in this mode, the BTES was not charged by heat gain from the cooling system. As a result, the heat extraction rate in heating mode in Load 2 was higher than Load 1. Energy balance of the low-temperature side of the system on a monthly scale in Figure 5 shows the magnitude of energy flows. For Load 1, in the first four months as well as the last three months, the evaporator load was higher than the cooling demand, because of the higher load on the condenser. In the rest of the year, higher cooling demand and surplus heat of condenser recovered the BTES energy levels to positive. For Load 2, in the first and the last third of the year when the system was operating in heating mode, the heat rate of the BTES was negative and in the second third of the year, cooling system injected energy to BTES in free cooling mode operation. Monthly share of the heating energy supplied to the heating system is shown in Figure 6. HWST was able to cover a maximum of 75% and a minimum of 45% at the monthly demand for Load 1. In warmer months due to lower accumulated energy in the HWST, the peak load contribution was higher. For Load 2, due to operation of the system for several consecutive hours in free cooling mode, the HWST was not charged enough to cover the demands of 6 th and 7 th month, while in the rest of the year the HWST covered between 55% to 92% of the monthly heat demand. Actual hourly loads on the evaporator and the condenser, as well as actual compressor power, are shown in Figure 7. Following the temperature reduction of BTES that was shown in Figure 4, as the temperature of BTES was decreased the heat pump performance was decreased with the same trend. Given the size of heat pumps chosen for the system, the heat pump operated in the full load for most of the year. The lowest partial operation hours were mostly in the heating mode. This was due to reducing of the BTES heat extraction capacity and high charging capacity of the water tank. For Load 2, in the heating mode operation, heat pump operated in the full load for the first part of the year. In the last 1000 hours of the year the performance of the heat pump was sharply decreased because of the low BTES temperature as can be seen in Figure 8.

Sensitivity analysis
The result of the sensitivity analysis is depicted in Figure 9. Among the evaluated parameters, only the impact of borehole pipe distances and the ground heat extraction rate on the volume of BTES were considerable.
The dependence of the effective length of the borehole pipes on the ground heat extraction rate was the only significant uncertainty for both Load 1 and Load 2. However, for Load 2 temperature of water at the condenser side and ∆ caused approximately 10% change in the total length of the borehole. The reason for this dependency was the change in evaporator operation regime that resulted in a change in the size of BTES.
The initial temperature of the BTES, the water temperature at the condenser side, demand magnitudes, and temperature difference at condenser side on the amount of heating energy covered by HWST were the most sensitive design parameters of HWST operated with Load 1. The HWST coverage was affected by more parameters for Load 2. Due to the utilization of a smaller HWST unit and therefore more prone to thermal disturbances thermal stratification, the volume of the HWST, the initial temperature of HWST, changed the heat coverage for Load 2 by nearly 10%. Borehole pipe distances and the brine temperature at the evaporator side had a larger impact on the HWST energy coverage for Load 2 comparing to Load 1. The reason for this change was the effect of these parameters on the BTES extraction rate and condenser actual load. Unlike Load 1, due to non-overlapping energy thermal demands for Load 2, the impact of change on the cooling demand and the water temperature at the condenser side was less significant.
Changes in both heating and cooling demands, the brine temperature at the evaporator side, initial temperature of BTES and temperature difference at the condenser side had the largest impact on , for both Load 1 and Load 2. For Load 2, borehole pipe distance, water tank stratification, and initial temperature of the water tank had a considerable impact on , .
The initial temperature of the BTES, magnitudes of demands and heat pump performance coefficients had the strongest impact on for Load 1. Regarding Load 2, the initial temperature of the BTES, the temperature at the secondary side of evaporator and condenser, borehole pipe distance and ∆ were the most impactful parameters.
The maximum heat extraction rate per meter from the BTES was affected by changing the initial temperature of the BTES, the brine temperature at the evaporator side, ground heat extraction rate, and cooling demand for Load 1. However, for Load 2 only the impact of the initial temperature of the BTES and the ground heat extraction rate were significant.

Long-term operation of the system
The simulation of the system operation in a 10-years period was performed on both Load 1 and Load 2 under the same design conditions for all the components. Besides the Load 1 and 2, the simulation was done for an alternative scenario with increased value by 0.20 for both Load 1 and Load 2.
As shown in Figure 10, for the base case the temperature of the BTES decreased for approximately six years before reaching a stable annual evolution for Load 1. This period was decreased for the alternative cases with higher to less than five years for Load 1 and almost three years for Load 2. Increasing the cooling to heating load ratio resulted in less reduction of ground temperature for both Load 1 and Load 2. The system performance coefficients are shown in Figure 11. The simulation with Load 1 resulted in a constant decrease in and , up to 20% and roughly 10%, respectively. The 10-year operation with Load 2 showed a 26% and 18% reduction of , however the rate of change after the sixth year with =0. 29 and fourth year with =0.49 was insignificant. Reduction in , for Load 2 was also evident. For =0.29, , was constantly decreased by up to 23% after the 10 th year, while with =0.49 the reduction nearly stopped after three years and 21% reduction.

Discussion
The model presented in this study aimed at the analysis of the sizing and thermal performance of the cooling and heating building energy systems. The focus was on the annual and 10-year study of the system with hourly time intervals. Therefore, some affecting aspects such as hydraulic flow distribution and control as well as the detailed structural design of components and distribution losses were not presented. The model was tested for two distinct combination and magnitudes of thermal loads. The results of the annual analysis reflected the expected thermal energy supply performance of the system within an acceptable range. Sensitivity analysis showed that the parameters affecting the sizing and performance of the system the most. Ground properties were found as an important factor in volume and effective length of boreholes. Temperature setpoints of the heat pump together with the intensity of thermal demands had a large impact on HWST coverage ratio. For Load 2 with lower total annual demands, the HWST energy coverage and , was more disturbed by the HWST design parameters such as thermal stratification and volume. Changing the initial temperature of the BTES largely affected indicating the dependence of system operation to the BTES temperatures. The temperature levels of the heat pump refrigerant had a small effect on , while the temperature set points at secondary sides had a larger impact. For Load 2 with a smaller BTES size, the effect of borehole distances was more significant due to an increased heat loss. The effect of changing cooling load, ground heat extraction rate, the initial temperature of the BTES, and the cooling brine temperature set point was considerable on the annual energy extracted from the BTES. Further investigation on 10-year operation showed that with a smaller cooling to heating ratio, the ground temperature will decrease more and consequently the and , will decrease on the long terms. Higher heat extracted from the cooling system together with lower load on the evaporator in heating mode would provide consequently higher available surplus heat and improved the temperature profile of BTES.

Conclusion
Sizing approaches for components of a heat pumpbased building energy system including heat pump, BTES and HWST as well as dynamic thermal energy balance of the system was introduced and implemented in MATLAB. The model was tested for two target users with completely different energy profiles and the sensitivity analysis on design parameters were applied. The result showed that the intensity of heating and cooling demands, as well as working temperature range of heat pump had a strong impact on system performance. Two definitions for the COP were introduced, one for heat pump stand-alone and one considering thermal energy storage units. The simulation of long-term operation showed a reduction of ground temperature and overall COP of heat pump as well as overall system efficiency. The ratio of annual cooling to heating load was defined to investigate the relation of load combinations with the long-term performance of the system. A stable and more efficient long-term operation was achieved for the system with higher . In general, efficiency reduction for the smaller system was sharper. However, on the long terms, the smaller system could reach the balance sooner than the larger system. As a further study, a detailed exergy analysis at the component level and considering the hydraulic aspect, as well as considering the integration of alternative energy sources such as solar thermal panels could give a clearer perspective of the performance efficiency of such systems.