A design tool for hybrid geothermal heat pump systems in heatingdominated buildings.
INTRODUCTION
Energy utilization in the built environment is of increasing concern, and geothermal heat pump (GHP) systems (also known as groundsource heat pump or Geoexchange[TM] systems) are now relatively well established as a means of significantly reducing energy consumption in space conditioning of buildings. This improvement in efficiency, however, generally comes at a higher first cost, which must be offset by lower operating and maintenance costs within an acceptable period of time to the building owner. As with most alternative energy systems, high capital cost is a significant barrier to market penetration. One of the main goals in the design a GHP system is to properly size the total length of the groundloop heat exchanger so that it provides fluid temperatures to the heat pump within design limits. Unlike with conventional heating and cooling systems, design of GHP systems requires some type of lifecycle simulation due to the thermal storage effects of the earth. Annual heating loads in a building are rarely balanced with annual cooling loads, and thus thermal responses of the ground throughout the building's life cycle must be considered. In heating dominated buildings, annual imbalances in the ground load will lead to progressively lower heat pump entering fluid temperatures, and in cooling dominated buildings, progressively higher heat pump entering fluid temperatures will occur. These excursions will result in the heat pump efficiency to progressively deteriorate and the equipment capacity to be ultimately compromised if the groundloop heat exchanger (GLHE) is not large enough. Borehole fields designed for buildings with relatively large annual ground load imbalances can often be excessively large and costly, making vertical closedloop GHP systems noncompetitive with conventional systems. The phenomenon of longterm temperature change in the subsurface due to GHP systems serving buildings with imbalanced annual thermal loads has given rise to the concept of the hybrid GHP system, where a supplemental component is utilized to effectively balance the annual ground loads. These systems permit the use of smaller, lowercost borehole fields, but their design adds to the complexity of the overall GHP design process because of the addition of another transient component of the system. For example, acceptable conditions for solar recharging of the ground in a heatingdominated building depends on solar availability and ground loop temperature. Consequently, hybrid GHP systems should be analyzed on an hourly basis in order to fully understand their behavior. The fundamental task in designing hybrid GHP systems lies in properly balancing the size of the supplemental component and the size of the groundloop heat exchanger, while optimizing the control of the supplemental component. Current engineering design manuals such as Caneta Research, (1995), Kavanaugh and Rafferty (1997), Kavanaugh (1998), and ASHRAE (2003), developed from research conducted since the 1980s, mention the potential use of hybrid cooling tower/fluid cooler GHP systems, but do not describe a design process for systems that utilize solar collectors for the thermal recharge of the ground. Proper and reliable design of hybrid GHP systems is quite difficult and cumbersome without the use of a system simulation approach. Further, without an automated optimization scheme coupled to the system simulation program, the design activity itself can become tediously impractical and time consuming. The use of system simulation for analyzing complex building systems is ever increasing, but the necessary computing resources are not at the disposal of every design practitioner, nor is their use always economically justified. As new technologies and design concepts emerge, design tools and methodologies must accompany them and be made usable for practitioners. Without reliable design tools, reluctance of practitioners to implement more complex systems can become a significant barrier. Therefore, the overall goal of the work presented in this paper is to develop a tool for the design of hybrid GHP systems in heatingdominated buildings that is useful to practitioners. The design tool is based on current "stateoftheart" system simulation methods, cast in a format that allows straightforward use. BACKGROUND AND LITERATURE REVIEW A review of the literature reveals a number of works that could be regarded as hybrid systems. Research and development regarding coupled solar and geothermal systems has largely been administered by the International Energy Agency (IEA) through working agreement'Energy Conservation through Energy Storage' (ECES), Annex 8 (Implementing Underground Thermal Energy Storage Systems) and Annex 12 (High Temperature Thermal Energy Aquifer and Duct Storage). In the United States, documented studies and reports in the literature dealing with hybrid GHP systems mainly pertain to cooling dominated buildings that utilize cooling towers. A summary of some of the heating dominated works is presented here. Recharging the ground with solar energy is not a new idea. Andrews et al. (1978) present qualitative evidence to suggest that coupling a ground loop to a solar collector and a storage tank could improve system performance and reduce initial costs. The need for earth heat transfer models was identified. Dalenback et al. (2000) present a preliminary study of a solarheated low temperature spaceheating system with seasonal storage in the earth for 90 singlefamily homes near Stockholm, Sweden. The total annual heating demand is 1080 MWh. The system performance was simulated using MINSUN and TRNSYS. The heating system is a subfloor heating system designed for 30[degrees]C supply temperature. Sixty percent of the heating load is to be met with the seasonal underground storage volume (with no heat pumps) and the remainder by electric heaters. The system consists of 3,000 [m.sup.2] of integrated roof solar collectors and a borehole storage volume of 60,000 [m.sup.3], whose temperature varies between 3045[degrees]C throughout the year. Seiwald and Hahne (2000) describe a solar district heating system with an underground storage volume located in southwest Germany. The solar contribution is about 50% of the total heat demand. At the time of their publication, the system consisted of 2,700 [m.sup.2] of flat plate solar collectors and a borehole storage volume of 20,000 [m.sup.3]. The underground storage volume is planned for expansion to 63,000 [m.sup.3] and to be capable of longterm heat storage at temperatures of up to 80[degrees]C. Regular operation of the system began in January 1999. Reuss and Mueller (2000) describe another solar district heating system with an underground storage volume in Germany. The total annual space heating and hot water demand is 487 MWh. The system consists of 800 [m.sup.2] of roofintegrated solar collectors, a storage volume of 9,350 [m.sup.3], and a buried 500 [m.sup.3] concrete water tank. Heat pumps are used to ensure the heating demand is met. The system performance was simulated with TRNSYS, and with government subsidies, the system operating cost was shown to be cheaper than conventional heating systems. Chiasson and Spitler (2001) present a system simulation approach to the design of a hydronic snow melting system for a bridge deck on an interstate highway in Oklahoma. A vertical, closedloop GHP system was designed to meet the heating requirement. The massive concrete bridge deck slab was proposed to be used as a solar collector in the summer months to thermally "recharge" the ground by circulating fluid from the bridge to the ground. An approximate 20% savings in the total GLHE length could be realized if the solar recharge design option is used. Chiasson and Yavuzturk (2003) examined the viability of using solar thermal collectors as a means to balance annual ground loads. A system simulation approach was used to model an actual school building with typical meteorological year weather data for 6 U.S. cities with varying climates and insolation. Conventional and hybrid solar GHP systems were sized for each case for 20 years of operation with a minimum design entering fluid temperature of 0[degrees]C. Lifecycle cost analyses indicated that the drilling costs must exceed the range of $19.68/m to $32.81/m for hybrid solar GHP systems to be economically viable. METHODOLOGY The design tool presented in this paper is based on correlation of dimensionless groups that are comprised of variables describing the design of geothermal heat pump systems. A total of 62 GHP system simulations were conducted in order to generate design data for subsequent correlation. This section describes the system simulations conducted and the development of the dimensionless groups. System Simulations The GHP system models were constructed in the TRNSYS modeling environment (SEL, 2000) using standard and nonstandard TRNSYS library components. The system configuration for the base cases (i.e. with no supplemental hybrid component) is shown in Figure 1, and Figure 2 shows the system configuration for the hybrid cases. In all simulations, the building loads were precalculated as described below and read from a file into a water source heat pump component model. [FIGURE 1 OMITTED] Component Model Description. Generic Buildings: Hourly heating and cooling loads on an annual basis for generic buildings were computed using eQuest (Hirsch, 2005), a graphical user interface coupled to the DOE2 simulation engine (York and Cappiello, 1981). Hourly loads were calculated for a school building and an office building in 10 cities in the United States and Canada: Albuquerque, NM; Anchorage, AK; Boise, ID; Dallas, TX; Halifax, NS; Houston, TX; Philadelphia, PA; Raleigh, NC; Phoenix, AZ; and Vancouver, BC. In addition, hourly loads were calculated for a church and an apartment building for Halifax, NS and Houston, TX. Detailed descriptions of the buildings can be found in Chiasson (2007). Groundloop Heat Exchanger Model: This component model is a modified version of that developed by Yavuzturk and Spitler (1999), which is a response factor (referred to as gfunctions) model based on the work of Eskilson (1987). Modifications used in this present study are described by Chiasson (2007), and include replacement of a steadystate borehole resistance with hourly modeling of thermal storage of the utube fluid, and modeling of transient heat transfer within the borehole grout. Heat Pump Model: This component model describes the performance of a watertoair heat pump that has been developed for GHP system simulations. Inputs to the heat pump model include the building loads, entering fluid temperature, and fluid mass flow rate. Quadratic curvefit equations to manufacturer's heat pump catalog data are employed to compute the heat of rejection in cooling mode, heat of absorption in heating mode, and the heat pump energy consumption. Outputs provided by the model include exiting fluid temperature, energy consumption, fluid mass flow rate, and coefficient of performance (COP). Flow Controls, Pump, and Heat Exchanger Models: The hybrid GHP systems were modeled as a primary/secondary loop system as shown in Figure 2. The primary circuit is the building loop plus ground loop heat exchanger and the secondary circuit is the hybrid component loop. The flow circuits are separated by a plate type heat exchanger (TRNSYS component type 5) modeled with a constant effectiveness. The system is modeled in this way to mimic design practices and to allow different fluids to be placed in each loop. For example, the solar collector fluid is typically an aqueous solution of 50% propylene glycol to avoid extreme freezing conditions. [FIGURE 2 OMITTED] Tee pieces, diverters, and pumps were modeled using TRNSYS standard library component models. Differential temperature control schemes are specified with the differential controller model, which was used to activate pumps and valves. A simple control scheme was used to simulate a variable frequency drive on the primary building loop pump; the flow rate for the current hour of the simulation is scaled to the peak flow rate according to the ratio of the current hourly load to the peak building load. A maximum pump turndown speed of 30% was assumed. Solar Collector Model: The solar collector component model is TRNSYS Type 73, a theoretical flat plate solar collector. A singleglazed solar collector was modeled, whose thermal performance is given by: q" = [F.sub.R][I.sub.t[theta]][([tau][alpha]).sub.[theta]]  [U.sub.L]([T.sub.ap]  [T.sub.a]) (1) where q" is the useful heat gained by the collector per unit area, [F.sub.R] is the heat removal factor, [I.sub.[tau][alpha]]is the total irradiation of collector at incidence angle [theta], [([tau][alpha]).sub.[theta]] is the transmittanceabsorptance product of plate at incidence angle [theta], [U.sub.L] is the heat loss coefficient, [T.sub.ap] is the temperature of absorber plate, [T.sub.a] is the ambient air temperature. The heat removal factor ([F.sub.R]) is further defined as: [F.sub.R] = m[c.sub.p]/[A.sub.c][U.sub.L][1  exp([A.sub.C][U.sub.L]F'/m[c.sub.p])] (2) where m is the fluid mass flow rate, [A.sub.c] is the collector area, [c.sub.p] is the fluid specific heat capacity, and F' is the collector efficiency factor. The incidence angle ([theta]) was computed at each hour of the year using solar angle correlations described in ASHRAE (2005). Convection losses are calculated hourly, based on wind speed from the weather file. The outlet fluid temperature from the collector is determined by an energy balance on the fluid. Control Strategies for Hybrid GHP Systems. Previous work by Yavuzturk and Spitler (2000) and Ramamoorthy et al. (2001) concluded that the most effective control strategy for hybrid GHPs was a differential temperature control, where the monitored temperatures were the exiting heat pump fluid temperature and the source/sink temperature used by the supplemental component. With regard to solar collectors, the source temperature is the temperature of the collector absorber plate. The control strategy operates the solar collector loop when the temperature differential between the collector absorber plate temperature and the heat pump entering fluid temperature exceeds a selected value (described below). The operation of the solar collector loop allows heat to be transferred to the ground loop for the thermal recharge of the ground storage volume. The differential control temperature selected was 5.00[degrees]C for the solar hybrid GHP system simulations. Chiasson (2007) conducted a sensitivity analysis of the control strategy on the annual energy transferred and on the solar collector pump run time. The sensitivity analysis included system simulations with double and half the differential control temperature relative to the base case. Varying the differential control temperature in the hybrid solar GHP systems resulted in no significant change in the solar energy collected or in the annual heat pump energy consumed. However, decreases in the solar loop run hours as a result of increasing the differential control temperature to 10.0[degrees]C were observed at 5.6% and 7.1% for the Boise and Vancouver School cases, respectively. These decreases in run hours directly impact circulating pump energy, and it would therefore be advantageous to utilize greater differential control temperatures than simulated here. GroundLoop Sizing, Base Case Simulations, and Parametric Runs. The base GHP systems were simulated using TRNSYS as depicted in Figure 1 for a 30year life with a time step of one hour, by iteratively adjusting the borehole depth until the minimum entering fluid temperatures (EFT) were within the target criteria. The maximum/minimum EFT to the heat pump is one of the critical design parameters in the sizing of any GHP system. According to International Standards Organization (ISO), capacities of geothermal heat pumps in closedloop systems are rated at 0.00[degrees]C for heating and 25.0[degrees]C for cooling. In heatingdominated climates, the minimum EFT is the more critical design parameter and for this study, 0.00[degrees]C was chosen as the minimum design heat pump EFT. Earth thermal properties and borehole construction details used in the base design cases and parametric runs are summarized in Table 1. The purpose of conducting the parametric runs was to obtain more data points for subsequent correlation. The base case simulates a borehole network in the earth with a thermal conductivity ([k.sub.ground]) of 2.42 W*[m.sup.1]*[K.sup.1]. The base borehole design is with a diameter of 114 mm with 25.4 mm nominal diameter Utube pipes evenly spaced within the borehole with standard bentonite grout (i.e. a thermal conductivity of 0.750 W*[m.sup.1]*[K.sup.1]). The thermal conductivity of thermallyenhanced grout was taken as 1.50 W*[m.sup.1]*[K.sup.1]. The volumetric heat capacity used in all simulations was 2.5 MJ*[m.sup.3]*[K.sup.1] and 3.9 MJ*[m.sup.3]*[K.sup.1] for the earth and grout, respectively. Table 1. Summary of Parameters Common to Base Case and Parametric Runs Case [k.sub.ground], Description Name W*[m.sup.1] *[K.sup.1] Base 2.42 114 mm diameter bore, 25.4 mm Utubes evenly spaced, standard grout k1 1.38 114 mm diameter bore, 25.4 mm Utubes evenly spaced, standard grout k2 2.08 114 mm diameter bore, 25.4 mm Utubes evenly spaced, standard grout k3 3.46 114 mm diameter bore, 25.4 mm Utubes evenly spaced, standard grout R1 2.42 114 mm diameter bore, 25.4 mm touching borehole wall, standard grout R2 2.42 114 mm diameter bore, 25.4 mm Utubes evenly spaced, thermallyenhanced grout R3 2.42 114 mm diameter bore, 25.4 mm Utubes touching borehole wall, thermallyenhanced grout Peak design flow rates were established at 0.16 L/s per peak refrigeration tons, and the borehole heat transfer fluid was specified as an aqueous solution of 20% propylene glycol. To simulate pumping power, a head loss of 30.5 m of water was assumed adequate, based on a design criterion of 0.0400 m/m of head loss per equivalent pipe length. The Minimization Function and Optimization of the Hybrid GHP Systems. The minimization (or objective) function is given in Figure 2, but has not yet been described. The goal of the supplemental component is to effectively balance the ground loads over each annual cycle. When this occurs, the annual minimum heat pump entering fluid temperature in heatingdominated buildings should be equivalent from yeartoyear, rather than progressively reaching the design temperature several years into the future. When the minimum temperatures are equivalent from year to year, the system is considered optimized, and the ground loop length is at a minimum. Thus, the objective function to be minimized is given by: Z = [summation][([T.sub.ann.min.]  [T.sub.design]).sup.2] (3) System simulations were executed for a 30year lifecycle. System optimization was conducted using GenOpt (LBNL, 2004), a generic optimization software package, employing the NelderMead optimization algorithm. Four cases were selected as candidates for the hybrid GHP system simulations. Solar hybrid GHP systems were designed for the school building in Anchorage, AK; Boise, ID; Halifax, NS; and Vancouver, BC. The TRNSYS input files were constructed as shown in Figure 2. The hybrid cases used all the same parameters as the base cases (i.e. not the parametric runs). For each hybrid GHP case, a solar collector array was simulated with a tilt angle equal to the site latitude, and the collector azimuth was due south. Simulations were run for 30 years with a timestep of one hour. The borehole depths and solar collector area were adjusted automatically by the GenOpt software package until the objective function (Equation 3) had been minimized using a convergence criterion of 0.001[degrees]C. The heat pump minimum design EFT used in the objective function was the minimum EFT observed from the equivalent base case. A constant heat exchanger effectiveness of 0.80 was used. The circulating pump in the solar loop was simulated by assuming a collector flow rate of 0.013 L*[s.sup.1]*[m.sup.2] of collector area after Beckman et al. (1977), and the heat transfer fluid was an aqueous solution of 50% propylene glycol. An example plot of the heat pump entering fluid temperatures for the base case and corresponding hybrid solar GHP case for the Boise, ID school is shown in Figure 3. A review of Figure 3 shows that the system has been optimized, as defined previously, since the minimum heat pump entering fluid temperatures are constant from year to year. [FIGURE 3 OMITTED] Development of Dimensionless Groups The Buckingham Pi Theorem (Buckingham, 1914) was used to identify dimensionless groups for GHP systems. The theorem states that if an equation involving k variables is dimensionally homogeneous, it can be reduced to a relationship among kr independent dimensionless products, where r is the minimum number of reference dimensions required to describe the variables. The dimensionless products are frequently referred to as "pi groups", given the symbol [PI]. The important variables in GHP design were identified as: * Total borehole length, L, m * Borehole spacing, H, m * Peak heating load, [q.sub.h], W * Peak cooling load, [q.sub.c], W * Seasonal heat pump COP, SCOP * Annual equivalent full load heating hours, [EFLH.sub.h] * Annual equivalent full load cooling hours, [EFLH.sub.c] * Annual load contribution from supplemental component, [Q.sub.s], J * Undisturbed earth temperature, [T.sub.g], [degrees]C * Design maximum or minimum entering heat pump fluid temperature, [T.sub.max] or [T.sub.min], [degrees]C * Thermal diffusivity of soil/rock, [alpha], [m.sup.2]*[s.sup.1] * Design life, t, s * Steadystate borehole thermal resistance per unit length of borehole, [R'.sub.b], m*K*[W.sup.1] Some of these variables were combined to simplify the problem. The peak building loads were converted to annual ground loads defined as: Some of these variables were combined to simplify the problem. The peak building loads were converted to annual ground loads defined as: [Q.sub.e] = 3600[q.sub.h][EFLH.sub.h]([[SCOP.sub.h]  1]/[[SCOP.sub.h]]) (4) for heating, and [Q.sub.r] = 3600[q.sub.c][EFLH.sub.c]([[SCOP.sub.c] + 1]/[[SCOP.sub.c]]) (5) for cooling, where [Q.sub.e] and [Q.sub.r] are the annual heat extraction and rejection loads (J), and the SCOP values were computed from the heat pump model. The undisturbed earth temperature was combined with the critical heat pump design entering fluid temperature to obtain a [DELTA]T variable. The critical design temperature is the minimum design entering heat pump fluid temperature in heatingdominated applications, and the maximum design entering heat pump fluid temperature in coolingdominated applications. The variables were simplified and expressed in dimensional form as: L, m [right arrow] L H, m [right arrow] L [q.sub.e, r], W [right arrow] [FLT.sup.1] [Q.sub.e], J [right arrow] FL [Q.sub.r], J [right arrow] FL t, s [right arrow] T [DELTA]T, [degrees]C [right arrow] [theta] [alpha], [m.sup.2] * [s.sup.1] [right arrow] [L.sup.2][T.sup.1] [R'.sub.b], m * K * [W.sup.1] [right arrow] [theta]T[F.sup.1] where [Q.sub.e,r] is the peak ground extraction or rejection load (i.e. the building load adjusted for the heat pump COP at maximum/minimum entering fluid temperature). The above list includes 9 variables and 4 dimensions, which would result in 5 dimensionless groups. However, it was considered desirable to limit the dimensionless groups to three in order to simplify the correlation and to facilitate visualization of the fitted surface. In the above list of variables, the design life and the borehole spacing were considered those that have the least design freedoms, and in the GHP system simulations described previously, the borehole spacing was kept at a constant 6.1 m (20 ft) and the design life was kept at a constant 30 years. This reduced the problem to 7 variables and 4 dimensions, resulting in 3 dimensionless groups. To describe heatingdominated cases, [Q.sub.r], [DELTA]T, [q.sub.e], and [alpha] were selected as the repeating variables, and the following [PI] groups were formed: ln([Q.sub.e]/[Q.sub.r]), ln([[R'.sub.b][([q.sub.e]).sup.[3/2]]]/[[([alpha][Q.sub.r]).sup.[1/2]][DELTA]T]), ln([L[([q.sub.e]).sup.[1/2]]]/[([alpha][Q.sub.r]).sup.[1/2]]) Each dimensionless group is expressed as a natural logarithm of the dimensionless group in order to keep the values of the dimensionless groups from getting excessively large, and to facilitate the visualization of the surface fitting procedure described in the next subsection. Surface Fitting of Dimensionless Group Data For each of the 62 simulation cases, [pi] groups were calculated from the simulation results, and an automated surface fitting software tool TableCurve 3D (SYSTAT Software, Inc., 2002) was used to fit an equation of a surface to the data set. The number of coefficients in fitted surface equations was limited to 10. The fitted equations were ranked in order of correlation coefficient. Prior to fitting a surface to the simulation data, some preprocessing of the data was necessary. A uniform grid was interpolated from the scattered simulation data using Surfer software (Golden Software, 2002) to avoid pitfalls of a surface equation being calculated by TableCurve 3D (SYSTAT Software, Inc., 2002) with unrealistic peaks and/or valleys in regions of the surface with sparse or lacking data. The method of kriging (Golden Software, 2002) was used to generate a uniform grid from the scattered data. The method of kriging is known as an exact interpolation method (as opposed to a smoothing interpolation scheme) that is mathematically similar to the method of least squares used in curve fitting (Golden Software, 2002). The kriging algorithm uses a semivariance model, which is a nonlinear weighting method based on the concept that sampled points that are closer to each other are more similar than those that are at greater distances. RESULTS AND DISCUSSION Correlating Equations of Dimensionless Groups Based on the simulation results, the dimensionless [PI] groups are defined as follows: [[PI].sub.1] = ln([Q.sub.e]/[Q.sub.r]) (6) [[PI].sub.2] = ln([0.98724[([chi]).sup.0.70345][R'.sub.b][([q.sub.e]).sup.[3/2]]]/[[([alpha][Q.sub.r]).sup.[1/2]][DELTA]T]) (7) and [[PI].sub.3] = ln([1.01252L[([q.sub.e]).sup.[1/2]][([lambda]).sup.0.49428]]/[([alpha][Q.sub.r]).sup.[1/2]]) (8) where [chi] is the ratio of the particular steadystate borehole thermal resistance (per unit length) to that of the base case (0.1823 m*K*[W.sup.1]) and, [lambda] is defined as the ratio of the particular ground thermal conductivity to that of the base case (2.42 W*[m.sup.1]*[K.sup.1]). The best fit surface equation found by the surface fitting software tool (SYSTAT Software, Inc., 2002) was a Taylor Series polynomial with a correlation coefficient of 0.998, and is expressed as: [[PI].sub.3] = a + b[[PI].sub.1] + c[[PI].sub.2.sup.[1]] + d[[PI].sub.1.sup.2] + e[[PI].sub.2.sup.2] + f[[PI].sub.1][[PI].sub.2.sup.[1]] + g[[PI].sub.1.sup.3] + h[[PI].sub.2.sup.[3]] + i[[PI].sub.1][[PI].sub.2.sup.[2]] + j[[PI].sub.1.sup.2][[PI].sub.2.sup.[1]] (9) where the constants a through j are a = 3.0857E+02 b = 9.2721E+01 c = 6.4197E+03 d = 6.7618E+00 e = 4.5793E+04 f = 1.3673E+03 g = 1.2561E01 h = 1.0935E+05 i = 4.9644E+03 j = 5.1576E+01 A contour plot of the surface described by Equation 9 is shown in Figure 4. Also shown in Figure 4, regarding the optimized hybrid GHP cases, are arrows depicting the path taken from the base case point to the corresponding optimal point of balanced annual ground loads. [FIGURE 4 OMITTED] A review of Figure 4 reveals that the hybrid GHP cases take similar paths, approximately perpendicular to contour lines, to their optimum point at [[PI].sub.1] = 0.00. The change in [[PI].sub.2] (i.e., [DELTA][[PI].sub.2]) can be related to the change in [[PI].sub.1] (i.e., [DELTA][[PI].sub.1]) to describe the path taken from a base point to the optimum point as shown in Figure 5. [FIGURE 5 OMITTED] Approximation of the Seasonal Heat Pump Coefficient of Performance The coefficient of performance (COP) of geothermal heat pumps is rated according to International Standards Organization (ISO) Standard 132561 at 0.00[degrees]C for heating and 25.0[degrees]C for cooling, and these data are readily available from manufacturer's catalogs. However, calculation of the annual ground loads for use in the design tool presented here requires knowledge of the seasonal heat pump coefficient of performance (SCOP), which is defined as the average COP of the heat pump over a season. Without a system simulation, this value can only be estimated because the COP is a function of the entering fluid temperature, which varies throughout the life cycle of the system. In order to approximate the SCOP, another set of correlations has been developed from the simulation results of this study. The ratio of the SCOP to the standard rated COP is correlated to the natural logarithm of the annual building loads ratio (i.e. not the grounds loads ratio, as this is not known without the SCOP) and the undisturbed earth temperature ([degrees]C). As would be expected, the system simulation results reveal lower heat pump SCOPs occurring in cases with more imbalanced loads at more extreme undisturbed earth temperatures. At first glance, this correlation appears oversimplified because two buildings, for example, could have equivalent annual loads ratios and undisturbed earth temperatures, but have different total ground loop lengths, and consequently have different maximum/minimum entering heat pump fluid temperatures. While this is true, the peak hour effects are apparently averaged out over the season within limits. Therefore, the limitations of the surface equation shown below are that ground loop lengths be sized to provide minimum entering fluid temperatures of 0.00[degrees]C[+ or ]3.00[degrees]C in heating mode and maximum entering fluid temperatures of 32.0[degrees]C[+ or ]5.00[degrees]C in cooling mode. These were temperature limits of the simulation data used to generate the contours. A surface was fitted to the simulation data using the same procedure as described to fit a surface to the [PI] groups. The equation of best fit for both the heating SCOP and cooling SCOP cases was a Taylor Series polynomial of the form: z = a + bx + cy + [dx.sup.2] + [ey.sup.2] + fxy + [gx.sup.3] + [hy.sup.3] + [ixy.sup.2] + [jx.sup.2]y (10) where z is the ratio of the SCOP to the rated COP, x is the natural logarithm of the annual building loads ratio, and y is the undisturbed earth temperature ([degrees]C). The coefficients of the surface equations are shown in Table 2. Table 2. Values of the Constants in the Fitted Surface Equations for Approximating Heat Pump Heating and Cooling Seasonal COP Heating SCOP Cooling SCOP a = 1.058E+00 a = 1.519E+00 b = 1.095E02 b = 2.393E02 c = 2.461E03 c = 1.155E02 d = 9.926E04 d = 9.887E04 e = 1.074E03 e = 1.174E03 f = 2.256E03 f = 3.369E03 g = 1.907E04 g = 3.238E04 h = 2.455E05 h = 2.863E05 i = 8.554E05 i = 1.387E04 j = 7.866E05 j = 1.089E04 To illustrate the use of Equation 10, consider a building where the annual heating load is four times greater than the annual cooling load, and the undisturbed ground temperature is 7[degrees]C. For a heat pump with a rated heating COP of 3.8 and a rated cooling COP of 5.5, Equation 10 predicts the ratio of heating SCOP to the rated heating COP (i.e., the z value) to be 1.08. Similarly, Equation 10 predicts the ratio of cooling SCOP to the rated cooling COP to be 1.46. This yields a predicted heating SCOP of 4.1 and a cooling SCOP of 8.0. If the annual heating load were two times greater than the annual cooling load (rather than four), Equation 10 predicts a heating SCOP of 4.2 and a cooling SCOP of 7.5. Analysis of Error The predicted values of [[PI].sub.3] using Equation 9 are plotted against the actual [[PI].sub.3] values in Figure 6. The average absolute percent difference between the predicted and actual [[PI].sub.3] values is 0.6%, with a maximum of 3.6%. Errors of similar magnitude are also observed in the calculation of heat pump heating and cooling SCOP using Equation 10. [FIGURE 6 OMITTED] A simplification made in the system simulations was to only deal with square borehole fields. Thus, there will be some error in adapting the results presented here to other borehole field patterns, because thermal interaction between the boreholes is, in part, a function of the placement of the boreholes relative to each other. To quantify this error, borehole fields were sized with varying configurations using GLHEPro software (Spitler, 1999) for the base cases of an extreme heatingdominated building (Anchorage, AK school) (Figure 7), and a building with wellbalanced annual ground loads (Halifax, NS office) (Figure 8). [FIGURE 7 OMITTED] [FIGURE 8 OMITTED] A review of Figures 7 and 8 shows that the greatest variation in total design ground loop length as a function of the borehole field pattern occurs for buildings with highly imbalanced annual ground loads. In changing from a 10x10 borehole field pattern to a 5x20 pattern, the same number of boreholes exist in each configuration, but a reduction in total ground loop length of 14% was observed for the Anchorage, AK school case due to the reduction in boreholetoborehole thermal interference. In changing from a 10x10 to an 8x8 borehole field pattern, the decrease in total ground loop length was not as great as in the 5x20 case, and was 6.8% for the Anchorage, AK toborehole thermal interference, but at the potential feasibility of increased drilling depths. For the balanced annual ground load case of the Halifax, NS office, the difference in total design ground loop lengths is less than 2% for the different borehole field patterns considered. This relatively low difference in design ground loop length is due to the relatively low boreholetoborehole thermal interference typical of balanced buildings; boreholetoborehole toborehole thermal interference is minimal because approximately equal quantities of thermal energy are extracted and rejected to the earth each year. Determining the Solar Collector Area The correlating equations described in the previous sections permit the calculation of total ground loop heat exchanger length, the quantity of annual solar energy required to balance the ground thermal loads, and the corresponding reduction of total ground loop heat exchanger length. To determine the solar collector area required to provide this annual heating load, the utilizability method is used. Details of this method are described in detail in Duffie and Beckman (1991). The concept is based on the premise that a critical solar radiation level can be defined below which no useful energy can be collected. Thus, utilizability is essentially a statistic which is a fraction of the total radiation that is received at an intensity higher than a critical level. The average radiation for the period can then be multiplied by this fraction to find the total utilizable energy. The method enables calculation of monthly quantity of energy delivered, given monthly values of incident solar radiation, ambient temperature, and heating load. The method of utilizability requires the inlet fluid temperature to the collector to be known. This temperature is unknown exactly without conducting some type of simulation, but can be approximated from the simulation results of this study by defining a monthly dimensionless inlet fluid temperature to the collector ([theta]): [theta] = [[[T.sub.collector in, avg]  [T.sub.earth]]/[[T.sub.earth]  [T.sub.min, design]]] (11) where [T.sub.collector in,avg] is the monthly average inlet temperature to the collector ([degrees]C), [T.sub.earth] is the undisturbed earth temperature ([degrees]C), and [T.sub.min,design] is the minimum design heat pump entering fluid temperature ([degrees]C). Thus, from Equation 11, a set of dimensionless typecurves can be generated (Figure 9) based on parameters that a designer would know (i.e. undisturbed earth temperature and minimum heat pump entering fluid temperature). Chiasson (2007) found that the calculated solar collector area by the utilizability method is relatively insensitive to the inlet collector temperature, within reason. For example, if the inlet temperature for each month was assumed to be the annual average inlet temperature, the increased error in collector area is only on the order of 1 to 2%. [FIGURE 9 OMITTED] In addition to the monthly average solar collector entering fluid temperature, the utilizability method requires design parameters that would be readily available to a practitioner. These include site latitude (degrees), collector tilt angle from horizontal (degrees), [F.sub.R][U.sub.L] (W*[m.sup.2]*[K.sup.1]) and [F.sub.R][([bar.[tau][alpha]]).sub.n]from collector test data, annual solar energy (J), monthly average daily solar radiation flux on horizontal (kWh*[m.sup.2]*[day.sup.1]), monthly average ambient air temperature. Note that [F.sub.R] is the heat removal factor, [U.sub.L] is the overall heat transfer coefficient for the collector (W*[m.sup.2]*[K.sup.1]), and [[bar.([tau][alpha]]).sub.n]is the monthly transmittanceabsorptance product at normal incidence. A comparison of the collector area calculated by the utilizability method to that calculated by the TRNSYS model revealed that the utilizability method tended to slightly overpredict the solar collector area. The smallest percent difference between the utilizability method and the TRNSYS model was 0.9%, observed for the Vancouver, BC school case, and the greatest was 6.6%, observed for the Anchorage, AK school case. The percent difference between the utilizability method and the TRNSYS model was 1.5% and 2.0% for the Boise, ID school and the Halifax, NS school, respectively. The above dimensionless temperature curves could only be considered valid to the situations modeled in this study. That is, with a heat exchanger effectiveness of 0.80 and the proportionate fluid mass flow rates and fluid properties. However, Beckman et al. (1977) provide a method to adjust the collector heat removal factor with other heat exchanger effectiveness values and fluid mass flow rates. For practical purposes, it can be assumed that the solar collector loop will have a smaller heat capacity rate than the ground loop in hybrid solar GHP systems giving the following expression: [[[F'.sub.R]]/[F.sub.R]] = [[1 + ([[F.sub.R][U.sub.L]]/[m[c.sub.p]])(1/[[epsilon].sub.c]  1)].sup.1] (12) where [F'.sub.R] is the solar collector heat exchanger efficiency factor [F.sub.R] is the heat removal factor, [U.sub.L] is the overall heat transfer coefficient for the collector (W*[m.sup.2]*[K.sup.1]), m is the collector fluid mass flow rate per unit area (kg*[s.sup.1]*[m.sup.2]), and [[epsilon].sub.c] is the constant heat exchanger effectiveness (). To apply the correction factor expressed in Equation 12 to the results of this study, [F'.sub.R]/[F.sub.R] would have to be first calculated for the values used in this study (i.e. a heat exchanger effectiveness of 0.80 and a volumetric flow rate of a 50% aqueous solution of propylene glycol of 0.013 L/s per [m.sup.2] of collector area) and then calculated for the particular values of interest. Then, the [F'.sub.R]/[F.sub.R] factor for use in the utilizability method calculations could be scaled accordingly. The Design Tool Algorithm An algorithm for a design tool for hybrid GHP systems using solar thermal collectors is shown graphically in Figures 10 and 11 and described below. The design algorithm requires input data that would be readily available to a designer. [FIGURE 10 OMITTED] [FIGURE 11 OMITTED] For standalone or hybrid GHP cases, input parameters required are a description of the building loads, the ground thermal properties, and heat pump details. Specific input data are shown in Figure 10. The first step in the algorithm is to calculate the heat pump seasonal COP from methods described above. The heat pump heating and cooling seasonal COP relative to its catalog rated COP has been correlated (Equation 10) to the average undisturbed ground temperature and the ratio annual heating to cooling loads of the building. The SCOP is then used to calculate the annual heating and cooling ground loads. The ratio of annual ground loads gives the value of the [[PI].sub.1] group as described by Equation 6. For positive [[PI].sub.1] values, the building is heatingdominated, and for negative [[PI].sub.1] values, the building is coolingdominated. The [[PI].sub.2] group is then readily calculated from Equation 7, and [[PI].sub.3] is then determined from the surface fit correlation (Equation 9), allowing calculation of the total ground loop length (L) from the [[PI].sub.3] expression described by Equation 8. The next set of calculations in the design algorithm (Figure 11) determines the required ground loop length for a hybrid GHP system, along with the area of the solar collector array. Additional required data are parameters that describe the heat exchanger and the solar collector array. Monthly average daily solar radiation on a horizontal surface, monthly ambient air temperature, and an estimate of monthly average collector inlet fluid temperature are also required. To calculate the new total ground loop length for a hybrid system, the annual energy contribution from the supplemental component is determined by setting [[PI].sub.1] equal to zero. The new value of the [[PI].sub.2] group is determined from the correlation shown in Figure 5. The calculation of the [[PI].sub.3] group and the total ground loop length is the same as for the standalone GHP case. Finally, the solar collector area is calculated by the utilizability method as described above. CONCLUSION An algorithm has been developed for a spreadsheetbased design tool for standalone and hybrid geothermal heat pump (GHP) systems that use supplemental solar thermal collectors in heatingdominated buildings. The algorithm has been developed from results of 62 detailed computer simulations, but such a design tool is not meant to replace system simulation for more detailed system analyses. Three dimensionless groups of key GHP design parameters were developed using the Buckingham Pi Theorem, and correlated using a fitted surface equation. The result is that from knowledge of the peak hourly heating and cooling loads, the annual equivalent full load heating and cooling hours, the thermal properties of the ground, the steadystate borehole resistance, the heat pump peak hourly and seasonal COP for heating and cooling, and the critical heat pump design entering fluid temperature, the total ground loop length can be calculated along with the quantity of annual energy required to balance the annual ground loads. With additional input parameters which would be readily available to designers, the area of a supplemental solar collector array in a heatingdominated building can be calculated, along with the reduced groundloop heat exchanger length. Solar collector array area is calculated using the utilizability method. Results from of a simplified design tool will expectably have some error when compared to an hourly system simulation due to simplifying assumptions and approximations made. The equations of best fit describing the correlating surfaces of the dimensionless groups result is a calculated total ground loop length error between 0 and 5%. This study assumed square borehole field patterns and an error analysis showed acceptable results in using the algorithm of this study to size the total ground loop length in closedloop, vertical borehole systems in other borehole field configurations. Solar collector area is calculated using the utilizability method as described by Duffie and Beckman (1991). This method provided a good match to the TRNSYS model results, with the smallest percent difference being 0.9%, observed for the Vancouver, BC school case, and the greatest being 6.6%, observed for the Anchorage, AK school case. NOMENCLATURE a,..., j = constants in surface fitting equations A = area, [m.sup.2] COP = coefficient of performance C = heat capacity rate, W*[K.sup.1] [c.sub.p] = specific heat at constant pressure, J*[kg.sup.1]*[K.sup.1] D = diameter, m EFT = entering fluid temperature F' = solar collector efficiency factor [F.sub.R] = solar collector heat removal efficiency factor [F'.sub.R]= solar collector heat exchanger efficiency factor H = borehole spacing, m k = thermal conductivity, W*[m.sup.1]*[K.sup.1] L = length, m m = mass flow rate, kg*[s.sup.1] q = heat transfer rate, W q' = heat transfer rate per unit length, W*[m.sup.1] Q = energy, J r = radius, m; R' = thermal resistance per unit length, m*K*[W.sup.1] SCOP = seasonal coefficient of performance t = time, s T = temperature, [degrees]C U = overall heat transfer coefficient, W*[m.sup.2]*[K.sup.1] Greek Symbols [alpha] = thermal diffusivity, [m.sup.2]*[s.sup.1] [CHI] = dimensionless steadystate borehole thermal resistance ratio [epsilon] = heat exchanger effectiveness [PI] = dimensionless pi group [theta] = dimensionless temperature [lambda] = dimensionless thermal conductivity ratio [rho] = density, kg*[m.sup.3] ([bar.[tau][alpha]]) = monthly average transmittanceabsorptance product Subscripts a = ambient air ap = absorber plate b = borehole wall c = cooling e = extraction g = undisturbed ground h = heating r = rejection s = supplemental component [theta] = solar incidence angle Overbar  = average condition REFERENCES Andrews, J.W., E.A. Kush, and P.D. Metz, 1978. A SolarAssisted Heat Pump System for CostEffective Space Heating and Cooling. Brookhaven National Laboratory Associated Universities, Inc. Sponsored by United States Department of Energy. ASHRAE, 2005. ASHRAE HandbookFundamentals, American Society of Heating, Refrigerating and AirConditioning Engineers, Inc., Atlanta, GA. ASHRAE, 2003. ASHRAE HandbookHVAC Applications, American Society of Heating, Refrigerating and AirConditioning Engineers, Inc., Atlanta, GA. Beckman, W.A., S.A. Klein, and J.A. Duffie, 1977. Solar Heating Design by the fchart Method. John Wiley & Sons, Inc., USA. Buckingham, E., 1914. On physically similar systems; Illustrations of the use of dimensional equations. Physical Review, Vol. 4, No. 4. Caneta Research, Inc., 1995. Commercial/Institutional GroundSource Heat Pump Engineering Manual. American Society of Heating, Refrigerating and AirConditioning Engineers, Inc., Atlanta, GA. Chiasson, A.D. 2007. Simulation and Design of Geothermal Heat Pump Systems. Doctoral Dissertation, Department of Civil & Architectural Engineering, University of Wyoming, Laramie, WY. Chiasson, A.D. and C. Yavuzturk, 2003. Assessment of the performance of hybrid geothermal heat pump systems with solar thermal collectors. ASHRAE Transactions, 109 (2). Chiasson, A.D. and J.D. Spitler, 2001. Modeling approach to design of a groundsource heat pump bridge deck heating system. Transportation Research Record, No. 1741, Advances in SnowRemoval and IceControl Technology: 207215. Dalenback, J.O, G. Hellstrom, S.E. Lundin, B. Nordell, and J. Dahm, 2000. Borehole heat storage for the anneberg solar heated residential district in Danderyd, Sweden. Terrastock 2000 Proceedings, 8th International Conference on Thermal Energy Storage: 201206. Duffie, J.A. and W.A. Beckman, 1991. Solar Engineering of Thermal Processes, 2nd Ed. John Wiley and Sons. Eskilson, P., 1987. Thermal Analysis of Heat Extraction Boreholes. Doctoral Thesis, University of Lund, Department of Mathematical Physics. Lund, Sweden. Golden Software, Inc., 2002. Surfer 8, Contouring and 3D Surface Mapping for Scientists and Engineers. Golden Software, Inc., Golden, Colorado, USA. Hirsch, J.J., 2005. eQUEST Energy Simulation Tool. Kavanaugh, S. P. and K. Rafferty, 1997. Ground Source Heat pumps: Design of Geothermal Systems for Commercial and Institutional Buildings. Atlanta: American Society of Heating, Refrigerating and AirConditioning Engineers, Inc., Atlanta, GA. Kavanaugh, S. P., 1998. A Design Method for Hybrid Ground Source Heat Pumps. ASHRAE Transactions 104(2): 691698. Lawrence Berkeley National Laboratory (LBNL). 2004. GenOpt 2.0, Generic Optimization Program, Version 2.0. Simulation Research Group, Building Technologies Department, Environmental Energy Technologies Division, Lawrence Berkeley National Laboratory. Ramamoorthy, M., H. Jin, A. Chiasson, and J.D. Spitler. 2001. Optimal sizing of hybrid groundsource heat pump systems that use a cooling pond as a supplemental heat rejectera system simulation approach. ASHRAE Transactions, 107(1). Ruess, M. and Mueller, J.P., 2000. Solar district heating with seasonal storage in Attenkirchen. Terrastock 2000 Proceedings, 8th International Conference on Thermal Energy Storage: 213220. Seiwald, H and E. Hahne, 2000. Underground seasonal heat storage for a solar heating system in Neckarsulm, Germany. Terrastock 2000 Proceedings, 8th International Conference on Thermal Energy Storage: 213220. SEL, 2000. TRNSYS Manual, a Transient Simulation Program, Version 15. Solar Engineering Laboratory. University of WisconsinMadison. Spitler, J.D. 1999. GLHEPro 3.0 for Windows Users Guide. School of Mechanical and Aerospace Engineering, Oklahoma State University. SYSTAT Software, Inc., 2002. TableCurve 3D, Automated Surface Fitting and Equation Discovery, Version 4.0. SYSTAT Software, Inc., Richmond, California, USA. Yavuzturk, C. and J.D. Spitler, 2000. Comparative study to investigate operating and control strategies for hybrid ground source heat pump systems using a short timestep simulation model. ASHRAE Transactions 106 (2):192209. Yavuzturk, C. and J.D. Spitler, 1999. A short time step response factor model for vertical ground loop heat exchangers. ASHRAE Transactions, 105(2): 475485. York, D.A. and Cappiello, C.C., 1981. DOE2 Engineers Manual (Version 2.1A), Lawrence Berkeley Laboratory and Los Alamos National Laboratory. DISCUSSION Jeffrey Spitler, Regents Professor, Oklahoma State University, Stillwater, OK: Are [[PI].sub.1] valves always expected to be 1? Andrew Chiasson: The value of [[PI].sub.1] is defined as the natural logarithm of the ratio of the annual energy extracted from the ground to the annual energy rejected to the ground. When the annual ground loads are equal, the value of [[PI].sub.1] is zero. [[PI].sub.1] is greater than zero for heatingdominated buildings and less than zero for coolingdominated buildings. Note that the definition makes use of the ground loads, which incorporates the building loads and the energy contribution by the supplemental component. While this sounds like a simple definition with little significance, it became quite significant when verified by the detailed TRNSYS simulations. That is, when the simulations were executed with the hybrid configuration shown in Figure 2, we found that [[PI].sub.1] travelled from its base case of imbalanced annual ground loads to a point of [[PI].sub.1] equal to zero with balanced annual ground loads when the GenOpt optimization routine was satisfied. Each hybrid case took a similar path toward [[PI].sub.1], which allowed their correlation. Richard A. Beier, Associate Professor, Oklahoma State University, Stillwater, OK: How did you determine the dominant dimensionless groups? With only three dimensionless groups in your model, I assume that not all of the parameters are represented in these three groups. How did you determine the dominant parameters and variables? Andrew Chiasson: All of the important design variables identified in the paper were incorporated into the dimensionless groups, except for the boreholetoborehole spacing, which is assumed constant at 20 ft (6 m). Some of the design variables were simplified. For example, the minimum (or maximum) heat pump design entering fluid temperature was combined with the undisturbed ground temperature to form a single DT variable. Also, as mentioned in the paper, square borehole fields were assumed and bestpractice heat pump flow rates were assumed. In forming dimensionless groups with the Buckingham Pi Theorem, the relevance of the groups can depend on selection of the repeating variable, and there is no unique set of groups. We systematically varied the repeating variables to find the most meaningful groups that resulted in an acceptable correlation. Better correlations are certainly possible, perhaps with four or five dimensionless groups. A.D. Chiasson, PhD, PE, PEng Associate Member ASHRAE C. Yavuzturk, PhD Member ASHRAE A.D. Chiasson is an assistant professor in the Department of Mechanical and Aerospace Engineering, University of Dayton, Dayton, OH. C. Yavuzturk is an assistant professor in the Department of Mechanical Engineering, University of Hartford, West Hartford, CT. 

Reader Opinion