{'fraction_non_alphanumeric': 0.03365751997221257, 'fraction_numerical': 0.013963181660298715, 'mean_word_length': 4.519746932515337, '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': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Analysis of short-to-long term heat flow in GSHP systems based on heat pump power

Abstract This paper presents a semi-analytical model based on the spectral element method for three-dimensional, short-to-long term heat flow in multiple borehole, multilayer ground source heat pump systems. The model is distinguished by its computational technique for expressing the input signal at the boundary of the borehole heat exchanger, giving rise to two important engineering features. First, the calculation can be conducted from seconds to years in a single run. This is achieved by discretizing the input signal at the inlet boundary of the borehole heat exchanger using a tailored fast Fourier transform with multiple time-stepping algorithm. Second, the calculation can be conducted using a Neumann boundary condition, instead of the commonly utilized Dirichlet boundary condition. This is achieved by mathematically relating the heat pump power to the heat flux at the inlet of the borehole heat exchanger, allowing direct use of the heat pump power signal as input instead of the inlet temperature. These features make the model computationally efficient that can readily be utilized for system design and included in inverse calculations. The two features are discussed in detail, verified against experimental measurements, and their functionality is highlighted by numerical examples.


Introduction
The ground source heat pump (GSHP) technology can play an important role in boosting the economy and improving the environment. It relies on energy gain from shallow depths, which are available nearly everywhere, and its operation produces minimal CO 2 emission into the atmosphere. As a result, this technology is booming, and geothermal engineers are endeavouring to make it more efficient and economic. Several European projects are specifically designed to develop more accurate and reliable computational models aiming to optimize the efficiency of this technology; see for example Horizon2010 Projects [1].
Presently, there are a considerable number of computational models describing heat flow in GSHP systems using different solution methods and adopting a wide range of physical and geometrical complexities. Cui et. al. [2] gave a comprehensive review on the available computational models for GSHP systems. Nevertheless, despite the presence of a large number of models and softwares, there are yet two important engineering features that need to be addressed in a more effective way: short-to-long term analysis of system performance, and analysis directly based on the heat pump power. This paper addresses these two features by introducing a computational technique enabling their implementation within the spectral element framework. Hereafter, we highlight the significance of these two features and how they are treated in this work.
Feature 1: Short-to-long term analysis is important to evaluate the system performance at different times of the days and years. Heating and cooling of buildings includes periods of switching ON and OFF, and during the GSHP lifetime, hourly, daily and seasonal periods of switching ON and OFF take place intensely. This operational requirement severely restricts the calculations such that if detailed transient analysis is sought, the long term performance of the system will be difficult to pursue, but if the global lifetime analysis is sought, the short term performance will be overlooked. As a consequence, the calculations are usually conducted for either short term or long term. This problem is in particular manifested in the numerical models and tools. COMSOL Multiphysics [3], FEFLOW finite element package [4], and, to a lesser extent, ABAQUS finite element package [5] and TOUGH2 finite difference code [6], are the most commonly utilized numerical tools for this purpose. These tools can provide advanced numerical facilities, but are severely restricted in conducting detailed short-to-long term analysis. The demand for a fine mesh due to the geometrical disproportionality of the system and the need for small time steps for the short term analysis make the CPU demands unrealistic for analysing the long term performance of the system, especially if the analysis is intended for daily engineering practice.
Analytical and semi-analytical models are not exemption and the majority of models are mainly suitable for either short term analysis or long term. Attempts to model short-to-long term performance of GSHP systems rely either on the Fast Fourier Transform (FFT) or the step response functions. Marcotte and Pasquier [7], and Zhang et al. [8] introduced, among others, semi-analytical models based on FFT. The FFT algorithm is rather efficient and employing it enables studying the performance of the system for as many time steps as required. However, the standard FFT is formulated for constant time stepping schemes. Depending on the size of the time step, the analysis can go for any desirable time span, yet, only one time step size is allowed: seconds, minutes, hours, days, months or years; not a combination of them. This characteristic restricts the standard FFT from utilization for short-tolong terms analysis; rather for either short term or long term.
Claesson and Javed [9], and Li et al. [10] introduced analytical models to calculate the temperature distribution in the system from minutes to decades. Claesson and Javed [9] adopted a two-step solution. For the short term (the first 100 h) they solve the conductive axialsymmetric heat equation analytically using Laplace transform; and for the long term, they use the line source solution. Li et al. [10], on other hand, developed a temperature response function that combines the composite-medium line-source solution and the conventional ILS and FLS solutions using the matched asymptotic expansion technique. More details on both models can be found in Li et al. [11]. Despite the effectiveness of these models, they are formulated based on the line heat source models which do not take the detailed conduction-convection heat flow in the borehole heat exchanger.
In this paper, we address detailed heat flow in an effectively 3D GSHP systems for short-to-long terms by tailoring the Fast Fourier Transform (FFT) algorithm to allow for the use of multiple time-stepping schemes. This feature facilitates the calculation of heat flow for any desired details, ranging from one second to years, done in a single run.
Feature 2: Analysis based directly on the heat pump power is physically more realistic than that based on a prescribed temperature at the inlet of pipe-in. This temperature is a result of the heat pump operation; i.e. not known a priori. In engineering practice, the heat pump power is prescribed as a heat flux in models solely solving the heat flow in infinite and semi-infinite domains, such as the line and cylindrical heat source models. For instance, the heat pump power can directly be prescribed as a constant heat flow rate per unit length, q, in the wellknown infinite line source model [12]: (1) in which T r t ( , ) is the temperature of the medium at radial distance r from the heat line source, s (W (m K)) is the solid domain thermal conductivity, (m s) 2 is thermal diffusivity and = 0.5772, Euler's constant.
Models which calculate heat flow in the borehole hole exchanger (BHE), the heat pump power cannot directly be prescribed. Instead, the temperature of the circulating fluid coming out of the heat pump and entering the inlet of pipe-in is prescribed. This temperature is calculated using the power equation: in which P (W) is the mean power of the heat pump, m (kg s) is the mass flow rate of the circulating fluid, c (J (kg K)) is its specific heat capacity, T (K) p in is the fluid temperature entering the heat pump and T (K) p out is the fluid temperature leaving the heat pump. T p in is considered equal to the temperature T out coming out of pipe-out of the BHE; and T p out is considered equal to the temperature T in entering pipe-in of the BHE, see Fig. 1. T in is utilized as the Dirichlet boundary condition to the inlet of pipe-in from which the temperature distribution in the GSHP system is calculated. The procedure to calculate T in from the heat pump power is given in Algorithm 1. This algorithm has been implemented in several numerical tools, see for example Al-Khoury et al. [13], Ozudogru et al. [14] and Rui et al. [15]. A similar algorithm is utilized in semi-analytical tools such as EED [16,17] and GLHEPRO [18,19]. Despite the practical use of this algorithm, it suffers the shortcomings of the step response function in terms of its computational inefficiency and the overlooking of the details within the time steps.

Algorithm 1 (Prescribing T in based on heat pump power).
In this paper we address the use of heat pump power as a Neumann boundary condition by formulating a mathematical equation relating the heat pump power to the heat flux at the inlet of pipe-in. This feature allows the direct use of heating and cooling design specifications to calculate the temperature distribution in the system. In this case, T in becomes unknown and needs to be computed, similar to T out and all other temperatures at any geometrical point in the system.

Theoretical background of the model
The above mentioned two engineering features, short-to-long term analysis and analysis based on the heat pump power, are implemented in the computer code STAND, which stands for Shallow geoThermal Analysis aNd Design. This code constitutes a semi-analytical spectral element model capable of computing detailed heat flow in effectively three-dimensional GSHPs constituting multiple BHEs embedded in multilayer systems, subjected to arbitrary temporal Neumann or Dirichlet boundary conditions.
The background theory of this model has been thoroughly presented in a string of publications, mainly Al-Khoury [20], BniLam and Al-Khoury [21,22] and BniLam et al. [23]. In Appendix A, we briefly present the governing equations and boundary conditions, together with the resulting spectral element equation necessary for utilization in the follow up sections.

Tailored fast Fourier transform
The fast Fourier transform (FFT) entails transforming a given function in its original time (space) domain to a frequency (wave number) domain and reversing it back. It computes the discrete Fourier transform (DFT) in an exceptionally efficient manner and saves significant CPU time and capacity.
For a given data function F t ( ) n with N samples, the DFT can be expressed as: in which t n is the discrete time sample, = k N 2 k is the discrete angular frequency and = k N 1, 2. .., 1. Inverse: Standard FFT requires N to be a power of 2 and the sampling rate (time step) constant for the whole signal period. The latter constraint is adequate in many applications in engineering, but for solving heat transfer problems in GSHP systems, it can be restrictive. In this kind of systems, the designer is interested in the detailed thermal response of the GSHP at certain periods of switching ON and OFF the system, while at the same time needs to know its life time performance. If a detailed analysis with a sampling rate of one second is conducted for the life time of the system, which is typically 20 years, it would require more than × 6.3 10 8 samples, which is practically unrealistic and can cause typical FFT numerical nuisances. If, on the other hand, an averaged analysis with a sampling rate of weeks or months is conducted, many important details on the hourly, daily and seasonal performance would be missing, making such analysis not representative of the real system. As such, the standard FFT can be useful to GSHP systems for either detailed analysis of relatively short periods, or averaged analysis of relatively long periods; not both at the same time.
The short time Fourier transform (STFT) is an elegant extension to the FFT, such that it enables dividing a time signal with a non-uniform frequency spectrum into segments, followed by computing the FFT for each segment separately [24]. The STFT, for segment j, can be expressed as: Inverse: Standard STFT requires equal interval segments and constant sampling rate to all segments. These requirements are adequate in communication and audio processing [24], but for GSHP systems, it is likewise restrictive.
To tackle this issue, we make use of the STFT idea to divide the time signal into segments, but tailor it to analyze segments with different time intervals and different sampling rates. Basically, we divide the signal into segments, use FFT with different time stepping rates (starting from small to large time steps), solve the system, and then link the reconstructed time domain of these segments via the Heaviside function, as illustrated hereafter. The tailored FFT is expressed as: where N j is the number of samples for segment j with time step t n j , and k j is the discrete angular sampling frequency of segment j, expressed as: Inverse: Using the Heaviside function, the reconstructed time domain of the signal is described as

N. BniLam and R. Al-Khoury
Applied Thermal Engineering 173 (2020) 114561 This equation expresses the time domain of a signal composed of N j j segments, put together. The algorithm of the tailored FFT can be summarized in Algorithm 2: Prescribe N j , t n j 3: Apply Forward tailored FFT, Eq. (7) 4: Conduct STAND calculations 5: Apply Inverse tailored FFT, Eq. (9) 6: Aggregate the results with j 1 segment, Eq. (10) 7: END DO The tailored FFT will be applied in the numerical example given in Section 6 to highlight the capability of the model to tackle the short-tolong time problems.

Heat flow analysis based on heat pump power
As indicated earlier, the heat pump power cannot directly be prescribed as a heat flux to models which simulate heat flow in the Utubes. Here, we derive an equation relating the heat pump power to the prescribed heat flux at the inlet of pipe-in, Eq. (A.7).
The power gained (lost) by the U-tube is expressed as Considering = m uA with A is the U-tube cross sectional area, Eq. (11) can be written as The Neumann boundary condition at the inlet of pipe-in, Eq. (A.7), is where (W (m K)) is the thermal conductivity of the circulating fluid. Substituting dT i of Eq. (12) into Eq. (13), and re-arranging, gives = q dz c u dP in (14) This equation entails that whether convective heat flux or conductive, the temperature at the inlet of the U-tube must be unique.
Integrating the left-hand side of this equation for a definite integral L 0 2 (with L being the BHE length, and L 2 denoting the length of Utube), and the right-hand side for P 0 , leads to which expresses the amount of heat transfer rate at the boundary of the U-tube when subjected to a heat pump power. The interplay between Eqs. (11) and (13) to produce Eq. (15) states that, in an adiabatic domain, the power gained by the U-tube (as a whole) is generated by a heat flux prescribed at the inlet of the U-tube. The heat flux at the boundary can be generated by any heat source, such as an electric heater or a heat pump.
As indicated above, Eq. (11) is valid for adiabatic processes, which occur without heat transfer between a thermodynamic system and its surrounding. Naturally, this is not the case for GSHP systems, as it is composed of multiple components and is in thermal contact with its surrounding soil mass. In mechanical engineering, the power of complex heat exchangers, such as those involving multiple tubes and shell passes, is calculated using a correction factor that implicitly incorporate the thermal interactions and losses among the heat exchanger components. Correction factors for several common configurations of mechanical heat exchangers are given in literature, see for example Pitts and Sissom [25]. Here, we apply a similar approach, by introducing a correction factor to Eq. (11), such that where is a correction factor reflecting the non-adiabatic process of the GSHP, and can serve to cover any anomalies in measured data and description of physics. This implies that the inlet heat flux, Eq. (15), must be modified to read: Normally, the heat pump power is not constant and varies with time due to electric voltage fluctuation and/or variation in its coefficient of performance (COP) [20], leading to: This time dependent heat flux signal is implemented in the spectral analysis by first transforming it to the frequency domain and the resulted signal is applied to q 1 on node 1 (Fig. A.2) on the right-hand side of Eq. (A.10).
To account for different BHE configurations, soil types and heat

N. BniLam and R. Al-Khoury
Applied Thermal Engineering 173 (2020) 114561 pumps, we determine the GSHP correction factor numerically by calibrating the thermal power using Algorithm 3.

Algorithm 3 (Tailored FFT & heat pump power prescription).
1: DO = j J 1 to (j: the time segment) 2: Prescribe N j , t n j 3: Transform the adiabatic form of the heat flux, Eq. (15), into the frequency domain and prescribe its magnitudes on node 1 in Eq. (A.10). 4: Conduct STAND calculations 5: Apply Inverse tailored FFT, Eq. (9) 6: Compute the power from Eq. (2) 7: Compare the specified heat pump power to the computed power, and calculate the correction factor, as Transform Eq. (17), into the frequency domain and prescribe its magnitudes on node 1 in Eq. (A.10). 10: Re-conduct the spectral analysis using STAND. 11: Apply the inverse tailored FFT, Eq. (9). 12: Aggregate the results with j 1 segment, Eq. (10) 13: END DO In Section 5 we validate this approach with experimental results and in Section 6 we utilize it in the numerical examples.

Model verification
Beier et al. [26] presented a well-instrumented and documented experimental test set-up for a prototype ground source heat pump. They introduced experimental results and provided digital data for heat and fluid flow in a single U-tube, representing a borehole heat exchanger, embedded in a horizontal 1.8 m × 1.8 m × 18 m sandbox. The borehole is 18 m long, centered along the length of the box, Fig. 2. They used an electric heater as a heat source. A fluid circulates through the closed loop and its temperature is measured at the inlet and outlet of the loop every 1 min. A flow meter is used to measure the fluid flow rate through the loop. Details of the experimental test set-up can be found in their paper, referenced above. Table 1 lists the involved sandbox and borehole parameters.
Two experiments were conducted: a continuous heat pump operation; and an interrupted heat pump operation, in which the heat pump is switched OFF after 9 h for 2 h, followed by switching ON till the end of the experiment. Both experiments ran for about 50 h. Using the provided digital data, we reproduced the experimental results by applying two calculation approaches: calculation based on prescribed T in , and calculation based on prescribed heat pump power.

Calculation based on prescribed T in
Prescribing the temperature at the inlet of pipe-in is a common practice in analysing heat flow in GSHP systems (see for example [27]). In this calculation approach we set the measured temperature at the inlet of pipe-in as a prescribed T in in the model, Eq. (A.8). Then, the temperature distributions along pipe-in, pipe-out and grout are calculated together with the temperature distribution in the surrounding soil mass. Fig. 3 shows the measured temperature at the inlet of pipe-in and the computed and measured temperatures at the outlet of pipe-out, for both experiments. The two plots in the figure show an excellent match between the experimental and computed results. Equally, Fig. 3b demonstrates the capability of the model to simulate accurately the switching OFF and ON the system in the interrupted operation experiment.
It can be noted that the model could also detect the drop of the inlet temperature in the interrupted operation experiment, occurred between 1 and 5 h, see the zoomed image in Fig. 3b. The model properly captured its effect on the output temperature, but it is not clear why this drop in temperature has not been reflected in the measured outlet temperature.

Calculation based on heat pump power
Calculation based on the heat pump power is conducted using Algorithm 3 (Section 4). Using this calculation approach the temperature distribution in the whole system is computed, including T in .

N. BniLam and R. Al-Khoury
Applied Thermal Engineering 173 (2020) 114561 anomalies in measured data and uncertainties in several parameters that resulted to good agreement between the experimental data and the computed temperatures at all stages: the transient at the beginning, the interrupted and the steady state.

Numerical examples
To demonstrate the model computational capability in simulating realistic cases, a numerical example illustrating heat flow in a typical ground source heat pump is introduced. A 2 × 2 BHE layout configuration constituting boreholes, 100 m in length, embedded in three soil layers exhibiting different thermal conductivities is simulated, Fig. 5. We studied three cases (Fig. 6 The material and geometrical parameters of the system are given in Table 2. The initial temperature in all components is assumed 12°C.
The heat pump works with varying power, depending on hourly, days and nights, and seasonal demands. Fig. 7a shows the power demand over the year, where it is assumed that this power represents the thermal energy gained from the heat pump, without the additional power usually gained from the compressor, neither the variations due to the COP of the heat pump.
The day starts at 6:00 (6 am) and ends at 22:00 (10 pm) and the night starts at 22:00 and ends at 6:00. Within these, the system is operating ON for 60 min and OFF for 20 min. Fig. 7b shows an example of the daily heat pump power for January. The year is divided into cold and warm months with varying thermal power demands. January till May, and October till December require heating; June requires neither heating nor cooling; and July till September are treated in two ways: with cooling (switching ON the heat pump in summer) and without cooling (switching OFF the heat pump in summer). The system is assumed to operate for 20 years. Any other operating scheme can also be considered.  The geometry is discretized using 3, 2-node spectral elements as schematically shown in Fig. A.2. The calculation is conducted using the tailored FFT with a multiple time stepping scheme, as the following: This scheme entails a very detailed analysis in the first day, a detailed analysis in the first year and a daily-averaged analysis for the rest. Other schemes can also be conducted, depending on the required details.
As indicated earlier, the model is capable of calculating the temperature distributions in all BHE components and in the soil mass in an effectively 3D space. Using the heat pump power option, even the temperature at the inlet of pipe-in, T in , is calculated. Fig. 8 shows the soil temperature distributions for the 3 cases: Case 1: One BHE is operating (Fig. 6) in two options: cooling ON in summer, and cooling OFF in summer. The temperature is computed at a 1 (0.1, 0.1), b 1 (2.5, 2.5) and c 1 (4, 4). The computed results for these three points are shown in Fig. 8a,  Fig. 8b, and Fig. 8c, respectively. Case 2: Four BHEs with × 5 m 5 m layout configuration (i.e. = 5 m, Figs. 5 and 6) in two options: cooling ON in summer, and cooling OFF in summer. The temperature is computed at a 2 (0.1, 0,1) and b 2 (2.5, 2.5). The computed results for these two points are shown in Fig. 8a and Fig. 8b, respectively. Case 3: Four BHEs with × 8m 8m layout configuration (i.e. d= 8m , Figs. 5 and 6) in two options: cooling ON in summer, and cooling OFF in summer. The temperature is computed at a 3 (0.1, 0,1) and c 3 (4, 4). The computed results for these two points are shown in Fig. 8a and Fig. 8b, respectively.

N. BniLam and R. Al-Khoury
Applied Thermal Engineering 173 (2020) 114561 soil temperature in the three cases is nearly the same for the cooling ON option, but for the cooling OFF, Cases 2 and 3 have clearly been affected. The soil temperature in these two cases is 3°C lower than if only one BHE is operating after 20 years of operation. At relatively far distances from the BHE, b 1 , b 2 , c 1 and c 3 , the effect of the borehole interactions are obvious. For Case 2, at point b 2 , the temperature difference after 20 years of operation between a single BHE and 4 BHEs is 4.3°C for the cooling ON option, and 7.0°C for the cooling OFF. For Case 3, at point c 3 , the temperatures were 3.3°C and 5.4°C, respectively. Fig. 8 also shows that the space between BHEs in a grid of neighbouring boreholes has a significant impact on the soil temperature. Likewise, cooling in summer is important for the soil mass to recover its thermal storage capacity. Fig. 9 shows the computed temperature variations for Case 2 at the inlet of pipe-in (T in ) and outlet of pipe-out (T out ) for two options: cooling ON in summer; and cooling OFF in summer. The figure shows the computed results for 20 years, together with the zoomed results for the first year and first day. These results demonstrate the computational capability of the model to calculate the detailed performance of the system for short periods, and its performance for long period, all done in a single run. Fig. 9 also shows an important physical observation. Operating the geothermal system for cooling in summer helps the ground to recover its heat and thus makes the system more sustainable. This observation has been discussed in detail by Zhao et al. [28]. In this specific case, the decline of temperature between year 1 and year 20 for the cooling ON is 1.8°C, but for the cooling OFF, it is 4.9°C. Fig. 10 shows the vertical temperature distributions for Case 2 in pipe-in, pipe-out, grout and soil at a 2 and b 2 in February and August of year 1, and those in year 20. The summer cooling is ON. The figure shows the model capability to simulate the detailed effects of the soil thermal conductivity on the thermal propagation and contraction at different months of the year. Fig. 11 shows the soil temperature snapshots at a × 10 m 10 m cross-sectional area surrounding the bottom of the boreholes It is worth mentioning here that although the model for a single borehole is axisymmetric, the use of the superposition principle has produced a significant spatial imbalance in the soil temperature distribution. However, by definition, the 2 × 2 BHE layout produces symmetric heat extraction rate for all BHEs. Any other layout such as 3 × 3 BHE or random distribution would produce unsymmetrical heat extraction rates. Kurevija et al. [29], Gultekin et al. [30] and You et al. [31] have discussed this issue in detail. This topic will be the focus of a coming publication.

Conclusions
This paper introduces a semi-analytical model for simulating heat flow in ground source heat pump systems, with emphasis on simulating short-to-long term heat flow using the heat pump power as input. The essential features of the model are: 1. It solves heat flow in effectively 3D geometries for any required temporal details. Detailed system operation from seconds to years can be handled in a single run. Hourly, daily and seasonally switching ON and OFF can be considered in the calculation. 2. It solves the system of equations using the heat pump power directly. Although in literature there are several models making use of the heat pump power to solve the system, they end up prescribing the Dirichlet boundary condition (T ) in at the inlet of the U-tube. Here, we prescribe the Neumann boundary condition (q ) in , derived from the heat pump power. The Neumann boundary condition allows using the heat pump power directly, but would be realistic if the COP and the compressor effect are taken into consideration. This subject will be the focus of a future work. Physically, the numerical example demonstrates that operating the heat pump in summer (i.e. cooling ON mode) helps the earth to recover its temperature, and hence increases the GSHP thermal efficiency.

Declaration of Competing Interest
There is no conflict of interest. in which the matrix on the left-hand side is the global stiffness matrix, with the superscript indices indicating the element (layer) number. The vector on the left-hand side is the degrees of freedom vector, indicating the nodal temperatures that need to be determined; and the vector on the right-hand side is the corresponding nodal heat fluxes. Readers interested in the formulation of the spectral elements for a multilayer GSHP systems are referred to BniLam et al. [23].  A2. A three-layer system and its spectral element discretization.