{'fraction_non_alphanumeric': 0.03730165687735601, 'fraction_numerical': 0.01564828614008942, 'mean_word_length': 4.4, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Hourly simulation of a Ground-Coupled Heat Pump system

In this paper, we present a MATLAB code for the hourly simulation of a whole Ground-Coupled Heat Pump (GCHP) system, based on the g-functions previously obtained by Zanchini and Lazzari. The code applies both to on-off heat pumps and to inverter-driven ones. It is employed to analyse the effects of the inverter and of the total length of the Borehole Heat Exchanger (BHE) field on the mean seasonal COP (SCOP) and on the mean seasonal EER (SEER) of a GCHP system designed for a residential house with 6 apartments in Bologna, North-Center Italy, with dominant heating loads. A BHE field with 3 in line boreholes is considered, with length of each BHE either 75 m or 105 m. The results show that the increase of the BHE length yields a SCOP enhancement of about 7%, while the SEER remains nearly unchanged. The replacement of the on-off heat pump by an inverter-driven one yields a SCOP enhancement of about 30% and a SEER enhancement of about 50%. The results demonstrate the importance of employing inverter-driven heat pumps for GCHP systems.

field at the BHE surface. The authors recommend to use the value given by that expression when averaged along the BHE length, which is called g-function. The g-functions are expressions of the time-dependent dimensionless temperature, averaged along the BHE length, due to a uniform and constant heat load which starts at the time instant t = 0. Simplified forms of the solution of Ref. [6], which has the form of a double integral, were proposed by Bandos et al. [7] and by Lamarche and Beauchamp [8]. By employing the results obtained in Ref. [7], Fossa [9,10] proposed simple approximate expressions for the g-functions based on the FLS model, which require low computational time and use empirical coefficients, determined through the analysis of different BHE fields.
Numerical simulations of BHE fields can be performed e.g. by means of the software Earth Energy Designer (EED), which is entirely dedicated to borehole heat exchangers and employs own g-functions, or through the software TRNSYS, which can perform an energy analysis of the whole building-plant system. However, as pointed out by Fossa and Minchio [11], TRNSYS simulations of BHE fields are not sufficiently accurate in the case of in line fields and/or unbalanced seasonal heat loads. Moreover, TRNSYS does not allow the simulation of heat pumps with inverter.
Recently, very accurate analytical expressions of the g-functions for the dynamic simulation of BHE fields have been obtained by Zanchini and Lazzari [12,13]. These g-functions, which are based on the Finite Cylindrical Source (FCS) model, are expressed in the form of polynomial functions of the logarithm of the dimensionless time. Through linear interpolations and superposition of effects, they can be applied for the dynamic simulation of any BHE field, provided that the ratio between length and diameter of each BHE is not lower than 500. In this paper, a code for the hourly simulation of GCHP systems, based on the g-functions obtained in Refs. [12,13], is presented. The code, executable through any programming language and here implemented in MATLAB, applies to monocompressor on-off and inverter-driven GCHP, used for building heating and/or cooling. Both the heat pump and the coupled BHE field can be simulated for several years. The code is employed to analyze the effects of the inverter and of the total length of the BHE field on the mean seasonal COP (SCOP) and on the mean seasonal EER (SEER) of a GCHP system designed for a residential house with 6 apartments in Bologna (North-Center Italy) with dominant heating loads.

Building and heat pump characterization
In order to characterize the building, the developed code requires as input data the hourly loads Eb during a whole year at the outlet of the generation subsystem (heat pump) for heating and for cooling. The building loads, which can derive from a dynamic simulation, must be set as negative for the heating season and as positive for the cooling season.
In order to characterize a mono-compressor on-off heat pump (ON-OFF HP), the required input data are the values, given by the manufacturer, of heat pump power, COP and EER in heating mode and in cooling mode, for a fixed value of the hot (or cold) water produced, Tw, and for different BHE fluid supply temperatures, Tfout. For inverter-driven heat pumps (IDHPs), these input values must be given in correspondence of the maximum, minimum and at least an additional intermediate inverter frequency.
Through interpolations of the manufacturer data, the heat pump power and COP in heating mode and the heat pump power and EER in cooling mode are expressed as second-order polynomial functions of Tfout, for each fixed value of Tw. In the case of inverter-driven heat pumps, a family of expressions for the heat pump power and a family of expressions for the heat pump COP (or EER) are obtained for each value of Tw, by varying the inverter frequency between the maximum and minimum value. The heat pump power and COP (or EER) expressions are stopped in correspondence of the cut-off temperature of the related operation mode. In fact, a minimum temperature of the BHE fluid could be fixed in winter (e.g. in order to prevent freezing, if just water is used as BHE fluid) and a maximum temperature in summer.

Borehole heat exchangers characterization
Regarding the borehole field, the BHE diameter (D), length (L) and thermal resistance per unit length (RBHE) must be set as inputs, as well as the number of boreholes (nBHEs) and their layout pattern and spacing.
The properties of the BHE heat carrier fluid (specific heat capacity at constant pressure, cpf, density, ρf, and volumetric flow rate, V̇f) and of the ground (thermal conductivity, kg, thermal diffusivity, αg, and undisturbed ground temperature, Tg) are also needed.
The mean dimensionless temperature at the interface between BHE field and ground is evaluated by the code by employing the analytical expressions of the g-functions obtained by Zanchini and Lazzari [12,13], which are given as functions of the logarithm of the dimensionless time, with dimensionless temperature, T * , and time, t * , defined as The term q0 in equation (1) is a reference thermal power per unit length. The mean value of the building thermal load during January (the month with the highest heating demand), divided by the total length of the boreholes, has been adopted for q0 in our case study.
The polynomial coefficients of the g-functions taken from Refs. [12,13] have been implemented in the code for several values of the BHE dimensionless length, L * =L/D, and for several values of the dimensionless radial distance from the BHE axis, r * =r/D.

Calculation of the GCHP system seasonal performance
In order to invoke the proper g-function coefficients, the code calculates the values of L * and r * required to evaluate the mean temperature at a BHE surface through the superposition of the effects in space. As example, for the case of three in line boreholes with L=105 m, D=0.15 m and a distance between adjacent BHEs of 6 m, the g-function coefficients are needed in correspondence of L * equal to 700 and r * equal to 0.5 (BHE -ground interface), 40 (central BHE -lateral BHE distance) and 80 (distance between lateral BHEs). If the g-function coefficients are not tabulated in correspondence of the required value of L * (or r * ), the closer lower value and the closer higher value of L * (or r * ) are considered by the code; then, linear interpolations of the results are employed.
At the beginning of the first hour of the simulated period, the BHE fluid is in thermal equilibrium with the surrounding ground; hence, the BHE fluid mean temperature (Tfm) and supply temperature (Tfout) are both initialized as equal to the undisturbed ground temperature (Tg).
For the i-th hour of the simulated period, the MATLAB code evaluates the season to which the hour belongs, reads the value of Tfout at the beginning of the i-th hour and calculates, through the heat pump second order polynomial functions, the heating (cooling) power that the heat pump is able to deliver and the corresponding COP (EER). For inverter driven heat pumps, a vector for the heat pump power and a vector for the corresponding COP or EER are obtained.
The energy delivered by the heat pump in the i-th hour, EHP (i), is equal to the product of the maximum heat pump capacity and the duration of one hour, thour (and is set as negative for the heating season), but if EHP (i) turns out higher than Eb (i), it is set equal to Eb (i).
For ON-OFF HPs, the values of the heat pump power and COP or EER in the i-th hour are known from the previous interpolations. The value of the heat pump power for IDHPs can be obtained dividing the module of EHP (i) by thour. If the result turns out lower than the heat pump power at the minimum inverter frequency, it is set equal to the minimum power (situation corresponding to on-off cycles). The corresponding COP, or EER, of the IDHP is then obtained by applying a second-order polynomial interpolation of the COP, or EER, vector as a function of the heat pump power vector.
The effective hourly values of heat pump COP or EER (COPeff (i) or EEReff (i)) take into account the heat pump efficiency decay in the case of on-off cycles, which occurs if, during the i-th hour, the heat pump power (at the minimum frequency for IDHPs) is too high compared to the building thermal needs. The values of COPeff or EEReff are evaluated according to the standards EN 14825 [14] and UNI/TS 11300-4 [15], multiplying the obtained COP or EER by the correction factor for on-off condition. The hourly value of the electric energy used by the heat pump is evaluated dividing EHP (i) by COPeff, or EEReff.
The thermal energy exchanged in the i-th hour between the borehole heat exchangers and the ground, Q (i), is evaluated as during the heating season, and as during the cooling season. Q (i) is negative if EHP (i) is negative, namely if heat is required by the building (extracted from the ground, during winter).
The mean value of the heat flux between BHE and ground per unit BHE length, q (i), is obtained as the ratio between Q (i) and the product of thour and the total length of the boreholes.
The dimensionless load amplitude of the i-th hour, A (i), is given by the ratio between q (i) and q0. At the end of the i-th hour, the dimensionless temperature T * m│r* (i), averaged along the BHE length, produced at the dimensionless distance r * from the BHE axis by a time-dependent dimensionless heat load, with steps of one hour and values given by the coefficients A, is: where the symbol g denotes the g-functions and t * hour is the dimensionless duration of one hour. By means of equation (5), T * m is calculated at r * =0.5 (BHE-ground interface), and at the dimensionless distances between the BHEs. If the g-function coefficients are not tabulated in correspondence of the required value of L * (or r * ), two hourly values of T * m are calculated by the code, in correspondence of the lower and higher closer available values of L * (or r * ). The actual hourly value of T * m for the given L * (or r * ) is then obtained through linear interpolation. The mean dimensionless temperature at the surface of a specific borehole of the field is evaluated as the sum of the value of T * m produced by the specific BHE at r * =0.5, and of those produced by the other BHEs of the field at their dimensionless distances from the specific BHE axis (superposition of the effects in space). The mean dimensionless temperature of the BHE field is evaluated as the arithmetic average of the T * m values of the different BHEs. The definition of dimensionless temperature of equation (1) yields the mean temperature of the BHE field at the end of the i-th hour, Tm,field (i). In the quasi-stationary approximation, the definition of BHE thermal resistance per unit length, RBHE, yields the BHE fluid mean temperature at the end of the i-th hour, Tfm (i): The corresponding BHE fluid supply temperature, Tfout (i), is obtained by subtracting to Tfm (i) half of the difference between the BHE fluid inlet temperature Tfin (i) and Tfout (i), namely The value Tfout (i) of the fluid supply temperature at the end of the i-th hour corresponds to that at the beginning of the subsequent ((i+1)-th) hour and is used by the MATLAB code to evaluate the heat pump performance at the subsequent hour, through a for cycle.
The hourly values of Tfm during the simulation period allow to check the long-term sustainability of the BHE field in the case of unbalanced building loads.
The seasonal performance of the ground-coupled heat pump during a selected year can be evaluated by means of the Seasonal Coefficient Of Performance (SCOP, ratio between the total energy supplied by the heat pump for heating and the corresponding electric energy used) and of the Seasonal Energy Efficiency Ratio (SEER, ratio between the total energy supplied by the heat pump for cooling and the corresponding electric energy used).
To check the long-term sustainability of a GCHP system, fast simulations of the whole system can be performed for several decades, by employing hourly simulations with the aid of auxiliary monthly simulations. In the case of monthly simulations, monthly building loads are given as inputs and a mathematical model very similar to that presented for the hourly simulation can be used.

Characteristics of the case study
A residential building which is undergoing an energy retrofit within the European project HERB (Holistic Energy-efficient Retrofitting of residential Buildings) has been selected to characterize the building heat loads. The house is located in Bologna (North-Center Italy), is composed of three floors with six apartments and has a total heated floor area of 282 m 2 . Figure 1 shows the power required by the building Pb at the outlet of the generation subsystem (heat pump), as a function of time, from October 1 st to September 30 th . The heating season has been set from October to April, included, while the cooling season from May to September, included. In figure 1 heating loads are considered negative, while cooling loads are considered positive, as required by the MATLAB code. Heating loads are dominant; the ratio between the annual thermal energy required by the building for heating and that for cooling is 2.32. The highest magnitudes of the power required by the building are 10.46 kW for heating and 9.32 kW for cooling.
The selected ground-coupled heat pump, used to provide heating and cooling to the building, is an inverter-driven brine-to-water unit. Water is delivered at 40 °C (return temperature 35 °C) during winter and at 7 °C (return temperature 12 °C) during summer. Data of heat pump power, COP and EER have been given by the manufacturer for BHE fluid supply temperatures (Tfout) from 5 °C to 18 °C in winter operation and from 15 °C to 35 °C in summer operation, at intervals of 1 °C, in correspondence of inverter frequencies equal to 110 Hz (maximum), 90 Hz, 70 Hz, 50 Hz and 30 Hz  The following ground properties have been set: Tg=14 °C, kg=1.8 W/(m K), αg=8.814 × 10 -7 m 2 /s. These properties correspond to a soil made of gravel and sand (55%) and clay (45%), partially saturated by groundwater, having a volumetric heat capacity equal to 2.042 × 10 6 J/(m 3 K).
The BHE field coupled to the heat pump is composed of three in line double-U boreholes, 6 m spaced, with diameter D=0.15 m and length L either 105 m or 75 m (corresponding dimensionless BHE lengths: L * =700 and L * =500, respectively). The BHE pipes, made of high-density polyethylene, have outer diameter 32 mm, internal diameter 26 mm and shank spacing 85 mm. The BHE is sealed with a grout having thermal conductivity 1.6 W/(m K).
The BHE fluid volumetric flow rate, V̇f, is 16 l/min. For the simulations with L=105 m, the BHE fluid is water, whose density ρf is 999.25 kg/m 3 and specific heat capacity cpf is 4.1896 kJ/(kg K) (properties at 14 °C). The corresponding winter cut-off temperature of the heat pump is 2 °C and the BHE thermal resistance per unit length (obtained through a numerical steady state simulation of the BHE cross section) is 0.0687 (m K)/W. Since simulations with L=75 m would cause too low water temperatures during winter, water is replaced by a mixture of water-glycol (monoethylenglycol 20%), whose ρf is 1032 kg/m 3 and cpf is 3.89 kJ/(kg K) (properties at 14 °C, from Ref. [16]). The corresponding winter cut-off temperature of the heat pump is -8 °C and the BHE thermal resistance per unit length is 0.0732 (m K)/W.

Results
Ten-year hourly simulations of the GCHP system have been performed with the selected IDHP, for the cases of L * =700 and L * =500. Figure 2 shows the hourly values of the BHE fluid mean temperature, Tfm, obtained during the 10 th year of simulation with L * =700. Each year starts on October 1 st and ends on September 30 th . The first ten points of the curves in figure 3 represent the minimum values of the mean temperature of the BHE fluid (Tfm,min), reached in each of the 10 years, for the cases of L * =700 and L * =500. The other points have been obtained by the aid of monthly simulations, as will be explained in the following. An unexpected peak of Tfm,min can be seen in figure 3 in correspondence of the second year of simulation. This is due to the fact that the first year starts with the heating season (heat extracted from ground), which, unlike the heating seasons of all the other years, is not preceded by a cooling season (heat injected in ground), that would enhance the fluid temperature. From the second year on, each heating season is preceded by a cooling season and Tfm,min has a regular decreasing trend.
Comparing the curves of figure 3, it is evident that, with L * =500, the BHE fluid reaches lower values of Tfm,min, compared to the case of L * =700, mainly due to higher heat loads per unit BHE length.
To analyze the influence of the heat pump inverter on the system seasonal performance during winter and during summer, the simulations are repeated forcing the heat pump to work at its maximum frequency, as a mono-compressor on-off heat pump. Table 1 shows the values of SCOP and of SEER of the ON-OFF HP and of the IDHP, for the 10 th year of operation, with both the considered BHE lengths. The performance values of Table 1 do not consider the energy used by the closed loop pump. From table 1 it can be noticed that the increase of the BHE length yields a SCOP enhancement of about 5% for the ON-OFF HP and of about 7% for the IDHP, while the SEER remains nearly unchanged. Indeed, since cooling loads are low, an increase of the BHE length does not yield a significant SEER improvement. The replacement of the mono-compressor on-off heat pump by an inverter-driven one yields a SCOP increase of about 30% for L * =500 and of about 33% for L * =700, together with a SEER enhancement of about 51% for L * =700 and of about 52% for L * =500. The most significant improvements on the seasonal performance are obtained by introducing the inverter, whose effect is greater on the summer efficiency (lower building loads, higher number of on-off cycles) than on the winter efficiency (higher building loads, lower number of on-off cycles). Figures 4 and 5 show, respectively, the hourly trend of COPeff and of EEReff, for the period November-January and June-August of the 10 th year, with and without inverter, for L * =700.  A simulation of the GCHP system has been performed for 50 years by means of monthly simulations. Moreover, hourly simulations of the heat pump from the 11 th year to the 50 th year have been performed by assuming that the hourly values of the BHE fluid mean and supply temperatures during a given year are equal to the corresponding values obtained at the 10 th year through hourly simulation, properly translated. The translation constant is the difference between the mean annual value of Tfm at the 10 th year and that at the given year, both obtained through monthly simulations.
The minimum annual values of Tfm so-obtained, from the 11 th to the 50 th year, are reported in figure 3 for the IDHP with L * =700 or L * =500. From figure 3 it can be noticed that the decrease in Tfm,min from the 2 nd to the 50 th year is less than 0.5 °C for L * =700 and less than 0.6 °C for L * =500. Therefore, in spite of the building unbalanced heat loads, the studied GCHP systems do not reveal long-term sustainability problems.
The accuracy of the assumption of employing monthly simulations for the translation constant has been tested by comparing the difference between the mean annual value of Tfm at the 5 th year and that at the 10 th year, obtained through monthly simulation, with the corresponding difference obtained though hourly simulation. The relative discrepancy between the two differences never exceeds 3.3% in the examined cases. Moreover, the accuracy of the assumption of employing a mean annual value as translation constant has been tested by analyzing the difference between the hourly values of Tfm at the 5 th year and those at the 10 th year, which is plotted in figure 6 for the case of the IDHP with L * =700. The plot of figure 6 shows that this difference is always very close to its mean annual value; the maximum deviation from the mean value is 0.01 °C. The seasonal performance coefficients obtained with the hourly simulations for the 50 th year are reported in table 2.

Conclusions
A code for the hourly simulation of Ground-Coupled Heat Pump (GCHP) systems has been developed. The code, which is implemented in MATLAB, employs the g-functions obtained by Zanchini and Lazzari [12,13] and applies to GCHPs with or without inverter, used for building heating and/or cooling. It allows fast hourly simulations for several years and, with the aid of auxiliary monthly simulations, even for several decades of a whole GCHP system, composed by the heat pump and the Borehole Heat Exchanger (BHE) field.
The code has been used to analyze the effects of the inverter and of the total length of the BHE field on the SCOP and SEER of a GCHP system designed for a residential house in Bologna with dominant heating loads. A BHE field with 3 in line boreholes has been considered, with length of each BHE either 75 m or 105 m. The results have shown that the increase of the BHE length yields a SCOP enhancement of about 7%, while the SEER almost does not change. Employing an inverter-driven heat pump instead of an on-off one can yield a SCOP increase of about 30% and a SEER enhancement of about 50%.
The obtained results demonstrate the importance of employing inverter-driven heat pumps for GCHP systems.