A Free Online Tool for the Design of Borehole Heat Exchanger Fields in Geothermal Heat Pump Applications at 10 and 25-Year Horizons

The proper design of borehole heat exchanger (BHE) fields used in ground-coupled heat pump (GCHP) applications assumes a key role in preserving the correct operability of the system. Each design method of BHE fields in GCHP applications requires a series of inputs. Among these inputs, there are the undisturbed ground temperature, the borehole and ground thermal properties, the choice of the BHE field configuration, i.e., the number and spacing of BHEs, and the pipe geometry. One of the most important inputs, sensibly affecting the design of BHE fields, is the building heat demand expressed in terms of a monthly heat load profile over a year and related peak loads in heating and cooling. By setting the desired Coefficient of Performance (COP) according to the geothermal carrier fluid temperature limits that cannot be overridden at the end of the considered time horizon, it is possible to correctly design the BHE field.

The ASHRAE method (also adopted by the Italian Technical Standard, UNI 11466), originally introduced by Kavanaugh and Rafferty [1], is one of the simplest algorithms for BHE field design. This method depends on the above-mentioned inputs and requires a reliable model for calculating the Temperature Penalty (T_p), which accounts for the thermal response of the borefield configuration using the g-function approach. In this method, the annual building heat load profile is condensed into three simplified heat pulses applied to the ground: 10 years of the yearly averaged heat load, 1 month of the most demanding monthly heat load, and 6 hours for the peak heat load of the most demanding month.

The method evaluates the ground’s thermal response as the time superposition of the solutions corresponding to the three aggregated heat pulses employing the Infinite Cylindrical Source (ICS) model [2]. However, since the ICS model does not account for 2D and 3D heat transfer effects in the ground, the corrective parameter T_p is introduced to improve long-term accuracy.

Several refinements to the original ASHRAE method have been developed. Fossa et al. [3] maintained the original ASHRAE formulation in their ASHRAE-T_p8 method but significantly enhanced the T_p calculation. Their approach introduced optimized constants tailored to different BHE field arrangements, derived from a broad analysis of g-functions by Eskilson [4] representing realistic ground thermal responses. The optimized constants are used in the expression for the T_p evaluation in the ASHRAE-T_p8 method.

All methods that deviate from the standard ASHRAE approach, such as those proposed by Ahmadfard and Bernier [5] and Fossa et al. [3], refine the evaluation of the overall borehole field length by typically assuming a plant lifespan of 10 years. Kurevija et al. [6] and other authors highlight the rising need for design methods and tools to correctly design and guarantee longer-term performance for at least 25 to 30 years. The 25-year horizon has been used in the comparison between coaxial, single and double U BHE configurations in terms of thermal and hydraulic performance reported in the recent study by Brown et al. [7]. Comparisons between a fully discretized 3D numerical BHE model, a resistance-capacitance (RC) model, a g-function model, and a hybrid model — against real-world monitoring data are provided by Heim et al. [8]. Three novel approaches to size geothermal borefields (energy piles) within a mixed-integer-linear-problem (MILP) framework, including a complex g-function-based model, have been proposed by Blanke et al. [9]. Their method proves to be particularly useful for large-scale energy system optimization, integrating accurate BHE sizing directly into MILP programming models. The work of Magdic et al. [10] considers the BHE and the rest of the system (heat pump, circulation pump, building) allowing to show that reducing borehole depth (while maintaining heat exchange) can yield more cost-effective savings over a 10-year operational period, without mentioning the results of the same analysis in terms of costs and BHE thermal performance for a more extended time horizon. The present authors were the first to extend the time horizon for the yearly averaged heat load from 10 to 25 years, covering a plant operability horizon more suitable and appealing for a commercial or residential unit. Thus, they propose the ASHRAE-T_{p 825} method that incorporates updated sets of optimized constants. Both the ASHRAE-T_{p 810} and ASHRAE-T_{p 825} methods have been extensively validated through a wide range of design test cases and related scenarios, including comparisons with commercial software [11, 12]. The method is now available as a free web-based design tool named BHEDesigner8, accessible at www.geosensingdesign.org, offering comprehensive BHE field design for both 10-year and 25-year planning horizons.

Comprehensive descriptions of the ASHRAE-T_{p 810} and ASHRAE-T_{p 825} methods, in particular, how the optimized constants have been obtained and how different shapes of the building heat load profile can affect the BHE field design for the 10-year and 25-year horizons, have been widely discussed in already published papers by the authors [11, 12].

This paper presents the key features of the online design tool and highlights the main factors distinguishing 10-year and 25-year designs. In particular, it discusses how the two different time horizon designs can be affected primarily by the yearly and peak heat loads, and variations in the estimated borehole thermal resistance (R_b) and ground thermal conductivity (k_gr). As evidenced in some reported cases, relatively small variations in the estimated k_gr and R_b can have a consistent impact across both horizons. In particular, the k_gr and R_b and their variations strongly impact the BHE field design determining an increase (or decrease) in the overall borehole length evaluated for each horizon design considered and influencing the level increment of L when evaluated in the 25th year instead of in the 10th year; the spread between the evaluated overall borehole lengths in the two design horizons tends to remain largely unchanged, with some exceptions (related to variations in k_gr). Recently, Tarrad [13] provided a correlation for the prediction of the borehole thermal resistance for the preliminary design of the U-tube ground heat exchanger. Ground thermal properties, borehole properties, ground pump parameters and economic considerations on the design and application of deep BHE have been discussed by Chen and Tomac [14].

The novelty and significance of this study lie in the implementation of the ASHRAE-T_{p 810} and ASHRAE-T_{p 825} methods within the first freely accessible web application designed to determine the required geometry of a BHE field serving a GCHP. Both the ASHRAE-T_{p 810} and ASHRAE-T_{p 825} methods, incorporated into BHEDesigner8, have demonstrated significantly improved precision compared to the original ASHRAE method while retaining its simplicity. In addition, the effectiveness of the T_p8 approach has been demonstrated by comparing its provided designs with those obtainable by established commercial software tools for BHE field design, such as EED and GLHEPRO [15].

In the present paper, the use of BHEDesigner8 is stressed once more in a series of test cases for both the above-mentioned design horizons. All calculations using the BHEDesigner8 have been performed directly through the free web application developed by the authors, and compared to those provided by the EED software. The case studies included in this analysis cover a wide range of BHE overall lengths and related depths and BHE field configurations, encompassing installations found worldwide.

The general definition of temperature penalty T_p is provided by Eq. (5):

$T_p(L)=\frac{\dot{Q}_y \Delta \Gamma}{L k_{g r}}$          (5)

where,

$\Delta \Gamma=\frac{\mathrm{g}\left(\tau_n\right)}{2 \pi}-G\left(\tau_n\right)$          (6)

∆Γ considers the error due to the use of the G temperature transfer function of the ICS solution for the long-term time $\tau_n$ and Eq. (6) compares it with the corresponding “true” Eskilson-type borefield g-function [4]. The time horizon (measured in seconds) $\tau_n$ relates to 10 or 25 years. The complete sets of equations and suitably tuned constants of the ASHRAE-T_{p 810} and ASHRAE-T_{p 825} methods are extensively detailed in papers authored by the present research group [3, 11, 12], to which the reader is referred for further details. All the equations of both the ASHRAE-T_{p 810} and ASHRAE-T_{p 825} methods are the same, along with the entire calculation structure, which iteratively converges to the L (and the related T_p) for each BHE configuration and set of inputs.

According to Fossa et al. [3], the ILS solution, along with its exponential integral function E₁, is used to evaluate the auxiliary parameter θ₈. This parameter is crucial for computing the temperature penalty T_p. To evaluate the θ₈ parameter, the spatial superposition of temperature solutions related to 8 BHEs arranged in a regular matrix configuration around the central one is applied. This superposition provides Eq. (7), where, B and $\mathrm{B} \sqrt{2}$ represent the center-to-center spacings between the central and the surrounding 8 boreholes.

$\theta_8=\frac{\dot{Q}_y\left[E_1\left(F o\left(\tau_n, B\right)\right)+E_1\left(F o\left(\tau_n, B \sqrt{2}\right)\right]\right.}{\pi k_{g r} \cdot L}$          (7)

The temperature penalty T_p8 is then expressed as:

$T_{p 8}=\frac{\theta_8\left(a N_4+b N_3+c N_2+d N_1\right)}{N_{t o t}}$          (8)

where,

N_tot is the total number of boreholes in the field.

N_j is the number of boreholes surrounded by j other boreholes [3, 11, 12].

The suitably tuned and optimized constants a, b, c and d are proper for the ASHRAE-T_{p 810} and ASHRAE-T_{p 825} methods. Different sets of constants in the range 0.03 < B/H < 0.125 (H is the depth of each borehole) are provided for different borehole configurations (rectangular and non-rectangular ones) and used according to the horizon design considered (the reader is addressed for further details to studies [3, 11, 12]).

It has to be highlighted that in the optimal search analysis aimed at the minimization of an objective function suitably defined in terms of T_p, a reference BHE depth (H_ref) of 100 m has been employed for evaluating each reference and “true” g-function. Following the approach suggested by Fossa et al. [3], for BHEs with different depths, both the ASHRAE-T_{p 810} and ASHRAE-T_{p 825} methods provide the corrected Fourier number Fo* expressed by Eq. (9) to be used for the calculation of the E₁ function in Eq. (7):

$F o^*=F o\left(\tau_n, R\right) \cdot\left(\frac{H_{r e f}}{H}\right)$          (9)

In a regular matrix arrangement comprised of 8 BHEs around the central one, the value of R can be either B or $\mathrm{B} \sqrt{2}$.

Whereas the main aim of the present section is to give further examples in the comparison of the designs for the 10 and 25-year horizon and arbitrary building heat load profiles could be assumed, a typical North European building structure related to a Swedish multi-family house built in a historical period between 1996-2005, located in Stockholm, has been considered for the present calculations. The building energy needs for heating and cooling and the hourly heat load profiles are obtained from EnergyPlus simulations reported in detail in Fossa et al. [3]. These simulations consider the local climate (as Typical Meteorological Years) and local building envelope characteristics (by the Tabula webtool database).

The monthly building heat loads resulting from EnergyPlus simulations are reported in detail in Figure 2 and are used as input for the BHEs field design tools; these heat loads are estimated by setting the internal temperature of the building to be maintained in the range 20-26℃ all year long. Moreover, based on the hourly EnergyPlus profiles, the peak loads in heating and cooling are estimated according to the approach proposed by Cullin and Spitler [22].

2.png

Figure 2. Monthly and hourly building heat loads profiles for a Swedish multi-family house from EnergyPlus

As shown by Figure 2, the Swedish case is characterized by an unbalanced heat load profile (i.e., heating mode largely prevails over the cooling mode) throughout the year; this unbalanced heat load profile acts as the base for further design cases. In fact, a second case of sizing with the ASHRAE-T_{p 810} and ASHRAE-T_{p 825} methods embedded in the BHEDesigner8 considers a BHE field serving 10 multi-family Swedish houses.

In this manner, the design process is addressed to involve both rectangular and non-rectangular borehole configurations and to highlight how the overall borehole length is highly dependent on the magnitude of the monthly heat loads considered. The same design cases of a BHE field serving 1 and 10 multi-family Swedish houses have been repeated for 3 different k_gr (1.5, 3, 6 W/mK) and R_b (0.05, 0.1, 0.2 mK/W) values to have an estimate of how the design is sensitive to variations in the estimated borehole and ground thermal properties.

All these design cases have been performed to guarantee that at the 10th or 25th year of operation both the monthly averaged COP_ave and the COP_peak are equal to 4.5 in the heating mode (winter season) and 4 in the cooling mode (summer season) of the GCHP. The expected minimum (or maximum) average fluid temperature T_f,ave(τ_n) corresponding to the assigned COP_peak is set to 0.5℃ in the heating mode (minimum value at the evaporator side of the GCHP) and 37.5℃ in the cooling mode (maximum value at the condenser side of the GCHP).

All these design cases have been performed with both the ASHRAE-T_{p 810} and ASHRAE-T_{p 825} methods embedded in the BHEDesigner8 and compared with the reference results (REF) provided by EED.

The intrinsic nature of Eq. (1) suggests considering also the important effect of the ratio $\frac{\dot{Q}_h}{\dot{Q}_y}$ between the most demanding hourly heat load $\left(\dot{Q}_h\right)$ and the yearly mean load $\left(\dot{Q}_y\right)$. Considering the same building heat load profile related to 10 multi-family Swedish houses (i.e., at the same $\dot{Q}_y$ and $\dot{Q}_m$), additional cases have been performed by setting different peak loads in order to vary the ratio $\frac{\dot{Q}_h}{\dot{Q}_y}$ (that is always greater than 1 by default). In particular, in some cases, the ratio $\frac{\dot{Q}_h}{\dot{Q}_y}$ has been doubled and halved of the original base case value.

Figure 3 presents the comparison for the single BHE depths H as computed by the reference tool EED (for both the 10- and 25-year time horizon) and the ASHRAE-T_{p 810} and ASHRAE-T_{p 825} methods. The two sizing approaches consider the same borefield configurations and the comparison in terms of H helps the visualization of all data. The ground volumetric heat capacity has been kept equal to 3 MJ/m³K. The ground thermal conductivity and borehole thermal resistance assume the above-mentioned values. The original values $\dot{Q}_y$ are 7.98 kW and 79.8 kW for 1 and 10 multi-family houses, respectively. The original value of the ratio $\frac{\dot{Q}_h}{\dot{Q}_y}$ is 4.21; some cases having $\dot{Q}_y$ of 79.8 kW have been run with doubled and halved $\frac{\dot{Q}_h}{\dot{Q}_y}$. The values assigned to the borehole spacing for these cases are 8, 8.5, 9.1, and 10.5 m (only for preserving a comparison between EED and T_p8 methods in the range of validity on B/H). The BHE fields considered in the validation are characterized by 4 × 4, 9 × 9, 10 × 12 rectangular, and 1 × 8 in-line (non-RCT) configurations.

3.png

Figure 3. Comparison between borehole depths as computed by EED, the ASHRAE-T_p 810 and ASHRAE-T_p 825

Figure 3 shows an excellent agreement with the EED method for the ASHRAE-T_{p 810} and a good agreement for the ASHRAE-T_{p 825}, demonstrating that the proposed one is a sufficiently accurate method to give a reliable design also for the 25-year horizon. As shown by Figure 3, the deviations between EED and ASHRAE-T_{p 825} are smaller for smaller borefields (for example, the one for a single multi-family house, i.e., the 4 × 4 rectangular and 1 × 8 in-line configurations); on the contrary, it becomes relevant for larger borefields.

Figure 4 shows the comprehensive comparison between the overall borehole length L as computed by the ASHRAE-T_{p 810} and ASHRAE-T_{p 825} methods for all the cases covering different $\dot{Q}_y, \frac{\dot{Q}_h}{\dot{Q}_y}$, k_gr and R_b. The borehole spacing and the ground volumetric heat capacity have been kept equal to 8 m and 3 MJ/m³K, respectively, for all the cases reported in Figures 4-10.

The impact of different k_gr and R_b on the design obtained with the ASHRAE-T_{p 810} and ASHRAE-T_{p 825} methods is separately detailed in Figures 5-10.

4.png

Figure 4. Comparison between overall borehole lengths as computed by the ASHRAE-T_p 810 and ASHRAE-T_p 825

As shown by Figure 4, the increase in the overall borehole length evaluated at the 25th year is included between 3% and 26% (the average increase, including all the cases, is 14%). The lower increments are for cases with $\dot{Q}_y$ quite close to 0 kW (7.98 kW), which is typical for balanced building heat load profiles along the year (the design cases involving the BHE field for a single multi-family house, i.e., the 4 × 4 rectangular and 1 × 8 in-line configurations); notice that when $\dot{Q}_y$ is quite close to 0 kW, consequently T_p is also close to zero according to its definition.

As shown by Figure 4, the deviations between the evaluated overall borehole lengths in the two design horizons become quite considerable when $\dot{Q}_y$ far from 0 kW (the design cases involving the BHE field serving 10 multi-family houses, i.e., the 11 × 13, 12 × 13, 14 × 14 rectangular configurations).

This would demonstrate that the design not only depends on the shape of the building heat load profile (balanced or unbalanced) and the related ratio between the yearly, monthly and hourly heat loads (at the ground side) but, more in general, it strongly depends on the magnitude of $\dot{Q}_y$ which assumes a primary key role in determining L and its increment when evaluated in the 25th year instead of in the 10th year. The $\dot{Q}_y$ far from 0 kW (for example, 79.8 kW) implies that it is comparable to $\dot{Q}_m$ because the heating or cooling loads prevail during the year (the cases with the $\dot{Q}_y$ of 79.8 kW have the $\dot{Q}_m$ of 289.8 kW).

5.png

Figure 5. Percentage deviations in overall borehole lengths between the ASHRAE-T_p 810 and ASHRAE-T_p 825 (1 unit)

6.png

Figure 6. Percentage deviations in overall borehole lengths between the ASHRAE-T_p 810 and ASHRAE-T_p 825 (10 units)

More consistent or attenuated deviations between L_25yrs and the L_10yrs, at the same $\dot{Q}_y$, are due to variations in k_gr, in R_b and in the $\frac{\dot{Q}_h}{\dot{Q}_y}$ ratio. Figures 5 and 6 represent the effects of k_gr and of R_b for different magnitudes of $\dot{Q}_y$, at the same $\frac{\dot{Q}_h}{\dot{Q}_y}$ ratio. For both considered $\dot{Q}_y$ values, increasing the ground thermal conductivity significantly reduces the deviations between L_25yrs and L_10yrs, as suggested by Eqs. (1)-(5). Variations of R_b, similarly to k_gr, considerably impact both the time horizon designs. Variations of R_b, differently to k_gr, tend not to modify the level of increment in L when evaluated in the 25th year instead of in the 10th year; the spread between the evaluated overall borehole lengths in the two design horizons tends to remain largely unchanged by varying R_b.

The figures also confirm that, as one could expect, $\dot{Q}_y$ far from 0 kW tends to increase and make more remarkable the percentage deviation between the L_25yrs and L_10yrs, provided the $\frac{\dot{Q}_h}{\dot{Q}_y}$ ratio is not particularly high. Doubling the $\frac{\dot{Q}_h}{\dot{Q}_y}$ ratio, at the same $\dot{Q}_y$, tends to decrease the percentage deviation by about 5% (compare Figures 6 and 7 with the corresponding k_gr), independently of the value assumed by B/H. This is particularly true and even potentiated by higher values of k_gr, even when $\dot{Q}_y$ is far from 0 kW, as confirmed by Figure 7.

7.png

Figure 7. Percentage deviations in overall borehole lengths between the ASHRAE-T_{p 810} and ASHRAE-T_{p 825} (different k_gr and $\frac{\dot{Q}_h}{\dot{Q}_y}$ at the same R_b and $\dot{Q}_y$)

Figure 8 collects and compares the absolute overall borehole lengths at 10 and 25-year horizons obtained by setting the three different k_gr values and the same R_b, $\frac{\dot{Q}_h}{\dot{Q}_y}$, and $\dot{Q}_y$ (10 multi-family houses). Figure 9 has the same meaning focusing on the cases for which the $\frac{\dot{Q}_h}{\dot{Q}_y}$ has been varied assuming the three different values. Figure 10 is related to the case of 1 multi-family house unit and reports the same information as Figure 8.

8.png

Figure 8. Overall borehole lengths at 10 and 25-year horizons (different k_gr at the same R_b, $\frac{\dot{Q}_h}{\dot{Q}_y}$, and $\dot{Q}_y$)

9.png

Figure 9. Overall borehole lengths at 10 and 25 years for different k_gr and $\frac{\dot{Q}_h}{\dot{Q}_y}$ at the same R_b, and $\dot{Q}_y$ (10 units)

10.png

Figure 10. Overall borehole lengths at 10 and 25 years for different k_gr at the same R_b, $\frac{\dot{Q}_h}{\dot{Q}_y}$, and $\dot{Q}_y$ (1 unit)

Figures 8-10 highlight once more how strong can be the impact of the ground thermal conductivity and how its correct estimation is crucial on the design of the BHE field. For this reason, the correct preliminary estimation via TRT assumes a key role. Overall, Figures 8-10 show how doubling the values of k_gr, at the same $\frac{\dot{Q}_h}{\dot{Q}_y}$ ratio, tends to decrease the absolute overall borehole lengths evaluated at the 10 and 25-year horizon designs and the related percentage deviation between the two horizon designs by about 5%.

This research is a part of the Italian National Program PNRR - Raise PNRR - ECS00000035 “RAISE (Robotics and AI for Socio economic Empowerment)” - SPOKE 3 “Environmental Caring and Protection Technologies, towards a Zero Emission Environment” - Research Activity (Task) 3.2 “Technologies for energy storage”, CUP D33C22000970006.