# A new approach to modeling ground heat exchangers in the initial phase of heat-flux build up.

INTRODUCTIONCurrently, ground-source heat pump (GSHP) systems are used extensively in the United States for space conditioning in commercial buildings. Their market penetration has increased progressively over the last decade, and presently they may be system of choice in many situations. The key design parameters of the ground heat exchanger (GHX) system are the number of boreholes, the depth and diameter of each borehole, borehole spacing, and fluid flow rates, etc., for attaining a target entering water temperature (EWT) after a specified number of years of operation on a design day.

The total mass of the fluid circulating in the GSHP systems, comprising the fluid mass in the GHX U-tubes and in the building heat pump circuit, is not entered as a parameter in the system models used in the design calculations. A key step in the design calculation is based on the assumption of heat exchange across the U-tube surfaces through the virtual presence of either a cylinder source (CS) or a line source (LS) of heat in the ground. The CS model assumes a steady heat flux across the boundary of the borehole of a known diameter into the surrounding medium. The LS model is based on the sustained release of the heat of constant strength (q'/pc) around an infinite line in the medium. Obviously, in these models, the thermal capacity of the mass of the aggregate fluid in the building and the GHX loops play no role.

In reality, the immediate response of the system to a step heat input is the gradual buildup of the heat flux across the bore boundaries and a slow rise in the EWT that is damped by the aggregate thermal mass of the fluid. The flux build up phase, characterized by increasing heat flux across the bore boundary is followed by the steady flux phase when the conditions are appropriate for applying the CS/ LS solution. This early response characteristic of the GHX is important in the design process as the peak EWT is the result of the long-term and the short- term thermal response of the ground to the building's load history. In the absence of an analytical model that captures the short-term thermal response to a reasonable degree, the researchers and the designers have used empirically summarized results of experimental investigations or numerical simulations. In this paper, the authors approach the problem intuitively with the fluid being modeled as a "virtual solid" and building heat being released as energy generated within the borehole. A classical analytical solution developed through Fourier transform has been adapted to fit the problem definition. The behavior of the average water temperature has been studied through results of the solution in dimensionless form.

SHORT TIME MODELING CHALLENGES AND PREVIOUS WORK

Historically, analytical models for the actual GHX, comprising the U-tube, grout, and borehole ensemble (Figure 1) are based on single-tube representations. It is argued that after some time the asymmetric heat penetration to the medium caused by the hotter fluid in the down pipe and the relatively cooler pipe in the adjacent leg can be ignored. Consequently the temperature distribution in the soil around the U-tube may resemble a thermal field developed by a single pipe of equivalent diameter. Thus the system of U-tube and the fluid inside is modeled as a single cylindrical core surrounded by the grout and the heat transferred to the fluid from the building is modeled as heat generated inside the core.

[FIGURE 1 OMITTED]

In the longer time frame, modeling issues for the GHX are the thermal interference of the adjacent boreholes and the heat flow across the boundary limits both in the lateral and axial directions of the GHX field. In the early stage of the flux buildup after a unit step heat input, the fluid temperatures have generally been obtained through numerical procedures (Yavazaturk et al. 1999; Young 2004; Xu and Spitler 2006). The key challenges for any reasonable analytical representation of the heat transfer in a GHX, in the current framework of model development, are as follows:

* Determining an equivalent diameter for the core to enable application of the analytical equations.

* Accounting for the thermal mass of the core and of the circulating fluid.

* Accounting for the impact of the grout.

* Ascertaining the impact of the thermal short circuit between the adjacent legs.

The Equivalent Diameter

Both the CS and LS analytical solutions are built around a single vertical element of heat source in the borehole. To meet the single-element criteria, the two legs of the U-tube are merged into a single pipe of an equivalent diameter, and the heat transfer is assumed to take place across the boundary of this pipe. Generally, this equivalent diameter [D.sub.eq] has been expressed as [D.sub.eq] = C * [D.sub.U-tube] where C is the correction factor determined by the geometry.

The determination of the equivalent diameter has been subject of considerable research (Claesson and Dunard 1983; Kavanaugh 1984; Bose 1984; Mei and Baxter 1986; Gu and O'Neal 1998; Sutton et al. 2002). One approach (Sutton et al. 2002) has been to find a C such that the actual thermal resistance of the borehole grout matches the value derived using the familiar steady state resistance formula for heat transfer across a thick-walled hollow cylinder, as given below:

[R.sub.p] = [ln([r.sub.0]/[r.sub.i])/2[pi][k.sub.p]] (1)

In an experimental study, Mei and Baxter (1986) observed large scatter in the values of C and concluded that a constant value of C cannot be applied with a "high degree of confidence." Hellstrom (1991) in his seminal work did not use equivalent diameter for deriving [R.sub.b] (the thermal resistance between the fluid and the borehole boundary). However various expressions of equivalent diameter are incorporated in several GHX design software.

Thermal Mass of the Core and the Circulating Fluid

In the LS solution, there is no discontinuity in the problem domain, and the thermal mass of the core is the same as that of the medium geometrically displaced by the pipe. In the CS model, the core is not part of the problem domain and hence not considered as constant heat flux is applied through the boundary of a cylindrical surface as a step change. Thermal capacity of the fluid is taken into account in numerical models (Muraya 1994; Shonder and Beck 1999; Rottmayer et al. 1997; Yavazaturk and Spitler 1999; Xu and Spitler 2006). Usually in the 2-D grids of the numerical domain, several cells of the grid represent the fluid mass. The GEMS2D code developed at Oklahoma State University has a provision for increasing the thermal capacity in the fluid cells by incorporating a fluid multiplication factor (Young 2004).

The heat carrier fluid exchanges heat simultaneously in the GHX loop and the heat pump loop of the building. In the building loop, the heat pumps reject the heat in the coaxial heat exchangers through which the carrier fluid flows at a preset rate. For a given cooling load (heat rejection rate), the temperature rise of the fluid leaving the building loop is inversely proportional to the flow rate. In the GHX loop, this heat is transferred to the ground. The net temperature response is damped by the thermal mass of the total fluid. The analytical models that have considered the fluid thermal capacity include the CS model (Model 3) used by Gehlin and Hellstrom (2003) and the Borehole Fluid Thermal Mass (BFTM) model (Young, 2004). Gehlin and Hellstrom used the large time solution for the heating of a semi-infinite body by a cylindrical probe with a given thermal mass (probe model). The solution, discussed in Carslaw and Jaeger (1959) is based on Blackwell's (1954) work. The heat capacity of the core was computed over the entire borehole weighted on the basis of the occupied areas of fluid and grout across a section. The borehole thermal resistance ([R.sub.b]) was an unknown parameter determined from the thermal response tests. Gehlin and Hellstrom (2003) used three sets of data to test four models and concluded that the probe model poorly agrees with other models (LS and numerical) and overestimates both the thermal parameters of [k.sub.s] and [R.sub.b]. Using the same formula, Hellstrom (1991) also worked out the upper time limit after which the fluid capacity may be neglected. Young (2004) has described a model where the buried electric cable analogy was used, and the temperature expression reported in Carslaw and Jaeger (1959) was used to compute the water temperature. The result was compared with a numerical model (GEMS2D), but the formula results did not compare well in many situations. This approach took into account the fluid mass outside the borehole and recommended that increasing the fluid mass outside the borehole may be considered to decrease the temperature spike due to a short heat pulse.

The Grout

As shown in Figure 1, to comply with the regulations to prevent aquifer contamination, it is common practice to pour a sealing material (grout) around the U-tube after it is placed in the borehole. The set grout mass also ensures good contact with open face of the borehole and the outer walls of the tubes thus aiding heat transfer. The most common material used as grout is 20%-30% solids bentonite clay that has sealing properties. In the unsteady state, if a single-step solution for the water temperature (as has been developed herein) is not available, the fluid temperature is calculated based on Equation 2 as follows:

[[theta].sub.f] = [[theta].sub.bb] + q'[R.sub.b] = [[theta].sub.bb] + q'([R.sub.g] + [R.sub.f] + [R.sub.p]) (2)

where [[theta].sub.bb] is the temperature excess at the boundary of the borehole, [R.sub.f] and [R.sub.p] are the steady state thermal resistances (neglecting the thermal mass of the fluid or the U-tube wall), and R is the steady state thermal resistance of the grout. The R values for the pipe and the fluid are given as:

[R.sub.p] = [ln([r.sub.0]/[r.sub.i]mp)/2[pi][k.sub.p]] (3)

and

[R.sub.f] = [1/2[pi][r.sub.i][h.sub.i]] (4)

Hellstrom (1991) has described three methods for the computation of the aggregate [R.sub.b] under steady state conditions that have been accepted as the most representative as of now.

Once the steady heat flux condition is reached, the [[theta].sub.bb] can be computed using either the CS or the LS solution. In the heat flux buildup phase, analytical models for the time varying grout resistance models [R.sub.g](t) or the [[theta].sub.bb] are not available. Gu and O'Neal (1995) have reported a general solution using an orthogonal expansion technique for heat diffusion in a composite cylinder under constant heat flux (CS) condition. But the problem domain does not cover a semi-infinite outer layer, and independent numerical computation of the eigenvalues leads to difficulties. Very recently, Lamarche and Beau-champ (2007) reported a solution under a constant flux condition that is able to accurately represent [theta] at the boundary of a constant cylinder source in composite media. When the same thermal properties are used for the grout and the medium, the results are identical to the CS solution.

The Short Circuit

The short circuit heat transfer that takes place because of the presence of two pipes containing fluid at different temperatures in close proximity to each other is a vexing problem to the modelers and during actual operation. In the numerical modeling with two pipes in the problem domain, this issue is not addressed as the two pipes are assumed to be at the same temperature. In single-core analytical modeling or two-pipe numerical modeling in 2-D this issue cannot be analyzed.

"VIRTUAL SOLID" (VS) MODEL

In the GSHP system, the heat rejected by the individual heat pumps to the carrier fluid circulating in the building loop reaches the ground-pipe interface through a combined advection and convection mechanism and finally diffuses to the medium by conduction dominantly in the radial direction. In the analytical models, the temperature excess (Q) at a reference point inside the borehole or at its boundary with the medium is expressed as a function of the applied heat flux, thermal properties of the fluid, grout and the medium and other geometric parameters. The generalized, nondimensionalized form of this expression is as follows:

[k[theta]/q'] = G(Fo, p1, p2,...) (5)

where [p.sub.1] and [p.sub.2] are the model parameters. The temperature excess [theta] represents the change in the temperature with respect to the initial condition, i.e., the undisturbed ground temperature or the far field temperature. In the LS solution, the G values are given by the following:

G(Fo) = [k[theta]/q'] = [1/4[pi]][[integral].sub.x.sup.[infinity]][[e.sup.-u]/u] du = [-1/4[pi]]Ei( - x) (6)

where x = l/(4Fo) and Ei = exponential integral. The CS solution is given as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (7)

In the above solutions, the thermal mass of the circulating fluid does not enter into the G value. In the stage of flux buildup, the thermal mass of the circulating fluid and the grout is significant in relation to the thermal mass of the ground in the immediate vicinity. Thus only a portion of the quantity of the heat released to the carrier fluid diffuses slowly in the ground in the periphery of the borehole, while the rest is absoreed in the fluid resulting in increased fluid temperature. The flux across the borehole builds gradually, and a steady flux may not be assumed till the elapsed time corresponds to the Fourier number of [Fo.sub.sf] that is dependent on the thermal diffusivity of the media and the grout as well as the borehole diameter.

In the present paper, the circulating fluid is modeled as a virtual solid. For any single borehole the specific fluid mass (per borehole total mass of the fluid in the system including the fluid in the U-tubes) is made up of the fluid inside the U-tube as well as a proportionate share of the fluid outside the U-tube i.e., in the building loop piping and the heat pumps. For a single borehole, the part of the system outside the borehole may be modeled as a well-insulated tank representing the building with the corresponding portion of the specific fluid mass in it. If the building is in the cooling mode, a prorated share of the heat transferred from the building to the fluid may be modeled as heat received by the fluid inside the tank from a heater coil. In the next step, the tank is dispensed with, and the heat is now generated inside the borehole. The fluid in the borehole is replaced by a solid with the same thermal mass as the specific fluid mass. The solid may be visualized as having the same density of the fluid, but having a higher specific heat (LpR times the specific heat of the fluid). This permits accounting of the aggregate thermal capacity of the circulating fluid inside the U-tube cubic space. As we are not interested in the fluid flow parameters per se, but in the heat transport aspects, the solid is assumed to have infinite heat conductivity, and a thermal contact resistance with the U-tube walls. In the final stage of the model development, the two legs of the U-tube are replaced by a single cylindrical rod having same cross-sectional area as that of the two tubes replaced. As shown in Figure 2, this is the core in which the heat is generated internally representing the heat transferred from the building to the fluid for eventual dissipation to the ground. The heat transferred to the fluid is now modeled as heat generated within the solid uniformly over the length of the solid cylinder in the core, similar to heat generated in the core of a fuel rod in a nuclear reactor or by resistive heating in a heating rod.

[FIGURE 2 OMITTED]

The key dimensionless parameter that characterizes the heat transfer in this system is the Biot number (Bi) given by

Bi = ([h.sub.eq]L/[k.sub.s]). (8)

The Biot number is often useful for the analysis of conductive-convective heat transfer between a solid surface and a fluid and is the ratio of the thermal resistance of the two media exchanging heat across a surface. While L is the characteristic length of the system, the [h.sub.eq] is the equivalent film conductance computed on the outer surface of the tubes of the U-tube (Figure 1) with the thermal resistance of the pipe wall section accounted for and is given as:

[h.sub.eq] = [[[[r.sub.o]/[r.sub.i][h.sub.i]] + [[r.sub.o]/[k.sub.p]]ln([r.sub.o]/[r.sub.i])].sup.-1] (9)

An identical relationship was used by Deerman and Kavanaugh (1991) in their analysis. In the VS model, several choices for analytical models exist as this problem of heat transfer is well studied (Carslaw and Jaeger 1959). Each analytical model has its limitations. Therefore, to fit the particular problem to the problem domain of the solution, simplifying assumptions are required. The solution used here was developed by Blackwell (1954). Carslaw and Jaeger (1959) have discussed the simpler large and small time expressions of the same.

This model, however, enables the average fluid temperature of the core fluid mass ([theta]*) to be determined directly. This imagery is limited to the heat transfer process between the ground and the fluid for the single-pipe representation in 2-D. All of the currently available analytical models also have identical constraints. However, it will be shown that the solution can represent the heat transfer across the U-tube geometry accurately with some limitations.

Analytical Solution

The analytical model is based on the solution of the unsteady state heat equation in the cylindrical coordinates given as:

[[[partial derivative].sup.2][theta]/[partial derivative][r.sup.2]] + [1[partial derivative][theta]/r[partial derivative]r] = [1/[alpha]][[partial derivative][theta]/[partial derivative][tau]] (10)

One of the boundary conditions for the medium represents the semi-infinite geometry, i.e.,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

This solution is only for the homogeneous medium extending up to the boundary of the VS cylinder. The virtual solid generating heat is in direct contact with the medium, but with a thermal contact conductance, [h.sub.eq]. The heat balances in the core volume and on the contact boundary surface produce the other boundary conditions:

[-k.sub.s][[partial derivative][[theta].sub.s]/[partial derivative]r] = [q'/2[pi]r*] - [r*[rho]*c*/2][[partial derivative][[theta].sub.s]/[partial derivative]t] at r = r*

and

[-k.sub.s][[partial derivative][[theta].sub.s]/[partial derivative]r] = [h.sub.eq]([[theta].sub.s] - [theta]*)at r = r*

The analytical solution (Blackwell 1954) is as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (11)

where

P = u[J.sub.0](u) + 2([[u.sup.2]/Bi] - [S/S*])[J.sub.1](u), O = u[Y.sub.0](u) + 2([[u.sup.2]/Bi] - [S/S*])[Y.sub.1](u)

Bi = (hD*/[k.sub.s]), and Fo = ([alpha][tau]/[r*.sup.2])

D*, the diameter of the core, is the characteristic length of the system. S* is the thermal capacity per unit length of the core while S is the thermal capacity of the medium of same volume, i.e., S/S* is the ratio of the heat capacity of an equivalent volume of the medium to that of the core.

The integral cannot be expressed as an algebraic function, but may be evaluated for various thermal capacity ratios and Biot numbers that depend on the fluid flow regime and geometry factors. The VS solution as described above is a single-step solution in the nonsteady flux state for the borehole heat transfer as it gives the average water temperature change in a single expression. This is applicable to a homogeneous medium and can be used directly in situations where the borehole cuttings are used to fill the hole after the insertion of the U-tube, as described in the ASHRAE manual (1997).

The time to reach the steady flux condition can be derived in straight forward manner from Equation 11 with the following:

Percent of energy absoreed in the core to energy released in the borehole

= [[([rho]*c*[theta]*)][pi][r*.sup.2]/[[integral].sub.0.sup.[tau]]q'*du = [[pi]/Fo](S*/S)G(Fo, Bi, S/S*) (12)

Example Reference Conditions

Numerical computations with actual "real-life" parameters are required for further observations on the behavior of this analytical solution. The geometry parameters and the soil thermal properties of the clay soil (Lake Agassiz sediment) used in all the subsequent calculations (Figures 3-9) pertain to an actual existing installation on the campus of the University of North Dakota in Grand Forks. These are presented in Table 1. The heat transfer fluid assumed in these calculations is water. Thermal properties of two other representative media are also listed in Table 1, and have been used in the calculations for comparing the G value responses of different media (Figure 4).

Table 1. Example Parameters Item Property Value Unit (SI) Value Unit (IP) Borehole Dia 114.3 mm 4.5 in. Heat 40 W/m 41.6 BTU/h*ft Input Rate Mediums Lake k 1.1 W/m.K 0.636 Btu/h*ft.[degrees]F Agassiz Clay [rho]c 2750 kJ/[m.sup.3] 41.0 Btu/[ft.sup.3] K [degrees]F Moist k 2.1 W/m.K 1.213 Btu/h*ft.[degrees]F Sandy Soil [rho]c 1756 kJ/[m.sup.3] 26.2 Btu/[ft.sup.3] K [degrees]F Typical k 3.0 W/m.k 1.733 Btu/h*ft.[degrees]F Rock [rho]c 2400 kJ/[m.sup.3] 35.8 Btu/h*[ft.sup.2] K [degrees] F Fluid h 1700 W/[m.sup.2] 299.4 Btu/[ft.sup.3] K [degrees] F Water [rho]c 4182 kJ/[m.sup.3] 62.4 Btu/[ft.sup.3] K [degrees] F k 0.6 W/m.K 0.347 Btu/h*ft.[degrees]F Pipe-1 ID 27.4 mm 1.08 in. in.IPS SDR-11 OD 33.4 mm 1.31 in. k 0.3895 W/m.K 0.225 Btu/h*ft.[degrees]F [rho]c 1900 kJ/[m.sup.3] 28.4 Btu/[ft.sup.3] K [degrees] F Tube Nil Tubes in "hug" mode (i.e., touching) spacing

The [h.sub.eq] computed using Equation 9 is 109.2 W/[m.sup.2]*K (19.2 Btu/h*[ft.sup.2]*[degrees]F), and the Biot number has been determined as 6.63 (basis explained later) for the clay soil. The equivalent diameter is [square root of (2)] times the tube OD i.e., 47.2 mm (1.85 in). The characteristic length used in the Biot number is twice the OD of the tubes of the U-tube assembly for the zero tube spacing configuration of the U-tube.

In Figure 3, the G values (dimensionless temperature) are shown at different loop fluid mass ratios (LpR) for the clay soil and for the tube geometry described in Table 1. The minimum value chosen for LpR is 1. Although, clearly some fluid outside of the U-tubes will always be present, the unit value of LpR is considered to facilitate analysis of one extreme. The reduced temperature response for the system with higher fluid mass is clearly evident. When compared between the highest LpR of 3.0 and the lowest LpR of 1.0, the G value for the system with the higher thermal capacity is lower by about 0.04 at around 1 hour. In the reference condition, this would imply a difference of about 1.5[degrees]C (2.6[degrees]F) in the average water temperatures of the two systems after imposing a heat pulse of 40 W/m (41.6 BTU/h * ft). The G value curves begin to converge after Log (Fo) = 1, i.e., Fo > 10 (corresponding to about 4 hours in real time for the given situation) showing the reduced impact of the thermal capacity of the core. In Figure 3, the G values as implemented in the ASHRAE design manual (1997) for the specific situation is also shown and it appears that the G values nearly converge with the VS solution for Fo > 60 (corresponding to about 23 hours in the reference condition).

[FIGURE 3 OMITTED]

From Figure 3, it may also be ascertained that if the thermal capacity of the circulating fluid is considered in the computation of the average water temperatures as is done in Equation 11, the design length of the GHX gets reduced as compared to the lengths obtained from conventional methods. In the design procedure, the design lengths are obtained on the basis of the target average water temperature set by the designer (ASHRAE 1997). For a given block heat load, the damping effect of the thermal capacity of the fluid results in the system responding with average water temperature lower than that obtained from the classical methods. The exact quantum of this reduction will depend on several factors i.e., the hours of duration of the peak load, the [k.sub.s] value, the loop ratio, etc. and may range from marginal to significant. Young (2004) studied the sensitivity of the GHX design lengths on core thermal capacity and several other factors for two reference buildings with highly different load profiles. For a one hour peak load dominated church building, he obtained a maximum reduction of 13.6% in the design length when the thermal capacity of the fluid in the U-tubes was considered. When the thermal capacity of the fluid outside the U-tubes (assumed to have same mass as the fluid within the U-tubes) is also considered, he observed an additional maximum reduction of 16.1%. Another advantage of this new approach will be improved accuracy in the building simulation programs for determining the annual energy consumption. This may be also used in the development of hybrid systems where a supplemental heat rejector is used in combination with a GHX.

Solution Behavior in Different Media

G plots for a loop ratio of two for three different types of soil described in Table 1 under identical conditions of step heat pulse are presented in Figure 4. The same plot also shows the CS solution for one of the mediums (the sandy soil) based on the steady flux at the borehole (with dimensionless [R.sub.b] computed at 0.24 using the approximate multipole formula) and Fo values adjusted for the change in the radius parameter from [r.sub.b] to r*. The adjustment is a mere change of scale required for plotting the results of the CS solution in the same Figure. In the CS solution, the Fo is based on the length scale of [r.sub.b] while the x-axis of the Figure 4 relates to Fo based on r* * It is evident that the thermal responses of different soil types are significantly different in the range of Fourier numbers of interest. The CS solution is essentially a steady flux solution and is not expected to yield accurate results at the early stage. But differences between the CS solution and the VS solution at a loop ratio of two (a very realistic condition) persist even for Fo > 100. The CS solution in Figure 4 is based on the steady flux at the borehole of 0.117m dia (4.5 in.) while the ASHRAE curve of Figure 3 is based on steady flux at a narrower diameter--the equivalent diameter of 0.055m (0.18 ft) resulting in superior convergence with the VS solution in Figure 3. The interesting aspect of Figure 4 is the way the heat capacities of the medium affect the temperature response of the fluid. In both the CS and the LS solutions, the heat capacity of the medium enters into solution only through the Fourier number. In the VS solution, the S/S* ratio is a significant parameter for evaluating the integral, and the shape of the G plots for different soils is significantly different because of this influence.

[FIGURE 4 OMITTED]

Comparison with the LS Solution

The large time expression for the VS solution as given by Carslaw and Jaeger (1959) is given as follows:

[k[theta]/q'] = [1/r[pi]][P + [4/Bi] + [1/2Fo]{P + 1 - [S/S*](P + [4/Bi])}] (13)

where P = ln(4Fo)-[gamma].

It is interesting to note that the first term of the above expression P/(4[pi]) corresponds to an approximation of the LS solution, while the second term (1/4[pi]) * (4/Bi) represents the difference between the core temperature ([theta]*) given by the VS solution as against the core boundary temperature given by the LS solution. The VS solution diverges from the LS solution as there are other terms that continue to remain important over a longer time period. The limiting value of the additional terms is zero for large Fo and thus the VS solution converges to the LS solution at large Fo values (P/Fo can be expressed as a power series of Fo, and its convergence behavior may be observed). Clearly, the solution passes the test of limiting value. The large time algebraic expression for G (Fo, Bi, S/ S*) also converges with the actual solution with the definite integral term given by Equation 11 for Fo > 100.

NUMERICAL SOLUTION

The analytical VS solution has been compared with the results from finite element models (FEM) in two dimensions for a single-pipe core and also for the U-tube assembly. A commercial FEM package (ANSYS) was used, and the key features of the FEM model used are:

* Thermal solid 2-D element PLANE77 to represent the soil /grout medium.

* Single layer of 2-D infinite solid element, INFIN110, at the periphery of the solid model to represent the semi-infinite boundary of the system. This element uses a mapping function to translate a semi-infinite boundary from a bounded problem domain.

* The core of each leg of the U-tube is represented by a thermal mass element MASS 71, its real constant set to the thermal capacity of the fluid in the pipe. The thermal capacity of the tube wall is ignored.

* The heat transfer from the tubes to the ground is captured through a number of convection elements, LINK34, with each element representing a 2.5[degrees] segment of the tube wall.

For further analysis and comparison of the results from FEM and the VS solution, the core equivalent for the U-tube geometry, described earlier is used as this enables direct application of Equation 11. The core is the solid rod with a diameter [square root of (2)] times the actual tube OD of the single tube of the U-tube assembly with identical thermal capacity of the filled U-tubes. This dimensional relationship preserves the area of cross section of the U-tube geometry. The thermal capacity is proportionally increased for loop ratios greater than one. The solid of the core has very high thermal conductivity, and no thermal gradient is sustained inside the e-tube except at the perimeter of contact with the medium. The contact conductance of the core is the [h.sub.eq] of the single tube as given by Equation 9. Heat may be generated within the core at a predetermined rate for unit length. The problem domain of the FEM with the core at the geometrical centre is shown in Figure 5, The ANSYS solution with appropriate real constants for the core at the center of the domain is completely congruent to the VS analytical solution, with matching parameters given in Table 2.

Table 2. Matching Parameters for ANSYS with e-Tube Core and the VS Analytical Solution ANSYS VS Solution 1 Thermal capacity of the core S/S* = 0.975 (clay/water) (Mass 71) set to the same for two Biot number based on the core tubes of the U-tube filled with diameter([square root of(2)]times water actual tube OD). 2 [h.sub.eq] = 109.2 W/[m.sup.2]*K Bi = 4.69 computed on the basis (19.2 Btu/h*f[t.sup.2]* of length scale equal to diameter [degrees]F) based on the OD of of the core, the [h.sub.eq], and real tube applied to the core the K for the clay soil. 3 Values for the thermal properties [rho]c of clay in Table 1 and of the clay soil used in the water used for the S/S * material properties description computation.

[FIGURE 5 OMITTED]

The Biot number for the geometry is first computed using core diameter as the characteristic length, as defined in Equation 8. The plot comparing the results of the FEM with a single-tube core and the analytical solution is not included in this presentation, but it is hardly surprising that both methods produce nearly identical results as they describe identical physical phenomena. This, however, serves the useful purpose of validating the FEM.

DISCUSSION

Tube Spacing Variation

When the single equivalent core is replaced by the actual U-tube with variable tube spacing in the FEMs (as shown in Figure 6), it is observed that the results agree with the VS model G values for a specific value of the Biot number. The matching value of the Biot number value, in turn, depends, for a given medium, only on the tube spacing ratio (defined as the distance between the walls of the two tubes divided by the tube outer radius). In Figure 7, the zero-spacing FEM result for the clay soil is compared with the VS solution for a Biot number of 6.63, and a close match of the two solutions is observed. This value of the Biot number corresponds to a length scale twice the OD of a single tube. The same exercises have been repeated for the other two mediums with three different tube sizes and have been found to be true in all cases. This is indeed a significant observation, as this enables direct computation of the fluid temperature in a real U-tube geometry without resorting to either an equivalent diameter or a steady flux assumption. Further, the incorporation of the thermal mass of the aggregate fluid in the model opens up other interesting possibilities.

[FIGURE 6 OMITTED]

[FIGURE 7 OMITTED]

In Figure 8, the FEM solutions for the U-tube geometry in the same clay medium are compared with the VS model solution with a low and a high tube-spacing ratio. At the lower SpR, the two solutions match closely. However at the higher SpR, the two curves intersect each other, and values diverge both at the lower and the higher end of Fo range. Evidently, the VS solution is not able to capture the physics of the process any more. In Figure 9, the Biot numbers and the corresponding tube-spacing ratios are plotted for the clay medium. It is seen that the solutions match closely. until a tube-spacing ratio of about 1.5, corresponding to a Biot number of 100, is reached. This is further validated with three different sizes of tubes DR 11 IPS pipes of nominal sizes of 0.75 and 1.5 in the same medium.

[FIGURE 8 OMITTED]

[FIGURE 9 OMITTED]

Several meaningful observations may be made on the Biot number values determined through the process of matching the FEM results with the results of the analytical solution T he extent of the increase of the Biot number with increasing tube-spacing ratio depends on the Biot number for the zero tube-spacing ratio ([Bi.sup.0] ), i.e., tubes touching each other. The influence of the Biot number on the temperature gradients in the convection-conduction heat transfer situation is well studied (Myers, 1998). It is observed that when [Bi.sup.0] is in the range of 2-3 as may be the case of high conductivity soils with the fluid in tureulent flow, the VS solution may remain valid for a tube spacing ratio of even up to 4, and the Biot number can double from the zero spacing value. On the other hand, in low conductivity soils with [Bi.sup.0] values in the range of 6-7, the Biot number may increase by approximately 12-15 times even with the increase of a tube-spacing ratio of 1.5. The VS solution results deviate from the FEM results with SpR values beyond this threshold value. The increase of the Biot number is equivalent to an increased convective heat transfer coefficient. Thus increasing the spacing between the tubes is equivalent to an increased fluid flow rate for obtaining higher Reynolds and Nusselt numbers. The VS model consequently would facilitate quantitative cost-benefit analysis for the usage of spacers, increased pumping energy to increase h, etc. A full table needs to be constructed showing the complete range of practically achievable [Bi.sup.0] values and corresponding Biot number with changing tube spacing ratio. A truncated version of the same is shown in Table 3.

Table 3. Tube Spacing Ration and Corresponding Biot Numbers Typical Biot Numbers Tube Spacing Low [k.sub.s] Medium High Comments Ratio [k.sub.s] [k.sub.s] 0 6.63 3.49 2.43 [Bi.sup.0] tubes touching and (typical values) 0.25 9 0.5 11 4 2.8 0.75 15 1 20 5 3.2 1.25 30 1.5 100 7 3.5 2 NA 8 4.0 NA: not applicable 2.5 NA 10 4.7 3 NA 15 5.0 4 NA NA 5.0~ ~Marginal increase

It is to be noted that Table 3 values are typical values obtained from the actual FEM behaviors for the example configuration and mediums described in Table 1. Further, Table 3 also points out the limitation of the VS model as applicable for U-tube geometry. In the medium to high k soils, it would possibly cover all practical tube spacing, though the same may not hold true for the low k soils.

Composite Media

The VS solution applies only to a homogeneous media, and it certainly restricts its usefulness to some extent. As discussed earlier, the large time algebraic expression of the VS model was used by Gehlin and Hellstrom (2003) for modeling the composite media of the borehole. They considered the entire borehole with the grout as the core generating the heat and transferred the entire thermal resistance consisting of the contributions from the fluid, the pipes, and the grout to the borehole boundary. This approach did not work well except for much longer times. Young (2004) used the buried cable model as reported in Carslaw and Jaeger (1959) and met with limited success. The buried cable model envisages a single thermal resistance sandwiched between two perfect conductors: the core and the sheath. Like Gehlin and Hellstrom (2003), Young also transferred the entire thermal resistance to the sandwiched layer and equated the thermal capacity of the sheath to the thermal capacity of the grout. As the results were not quite agreeing with the results from the numerical model (GEMS2D), Young adjusted the model by allocating a portion of the thermal capacity of the grout to the core. This improved the fit, and this factor was termed as the grout allocation factor.

It is foreseen that, at the initial stages of a step heat input in the borehole, the thermal front is within the grout envelope and has not seen the borehole boundary. The Biot number applicable in the VS solution would, consequently, be determined by the tube-spacing ratio and [k.sub.g] and is different from the Biot number in the late stage. If[k.sub.g] < [k.sub.s], the early-stage Biot number may increase proportionally by the conductivity ratio ([k.sub.s]/[k.sub.g]) for the zero tube-spacing configuration. However, as seen in Table 3, the Biot number may increase quite disproportionately for an arm's length configuration (tube spacing ratio [much greater than]0). The change in the Biot number depends on the applicable values in Table 3 of the [Bi.sup.0] values for the medium and the grout. Thus the early-stage G values do not depend on either the borehole diameter or on the medium conductivity. Indeed, this is clearly observed in the G plots from the FEMs for the grouted boreholes (Figures 10 and 11).

In Figure 10, the G plot for a typical grouted borehole is shown with the material properties modeled after the example given by Young (2004). The borehole and the U-tubes are identical to the ones described in Table 1 except that the pipe spacing ratio is 0.19. The [k.sub.s] value is 2.5 W/m*K (1.44 Btu/ h*ft*[degrees]F), and the pc value is set at 2500 kJ/[m.sup.3]*K(37.3 Btu/[ft.sup.3]*[degrees]F).The grout chosen is standard 20% bentonite with a [k.sub.g] of 0.75 W/m * K * (0.433 Btu/h*ft*[degrees]F) and a [rho]c value of 3900 kJ/[m.sup.3]*K * (58.2 Btu/[ft.sup.3]*[degrees]F). In the first step, G values and the corresponding [Fo.sub.g] values are computed from Equation 11 with an S/S* value and the Biot number as applicable for the specific grout/tube- geometry combination. Since in the G plot of Figure 10, the y -axis values are computed on the basis of [k.sub.s], and the x -axis values on the basis of Fo corresponding to [[alpha].sub.s], the first-step G values and the [Fo.sub.g] are adjusted for plotting on G([k.sub.s][theta]/q') vs.[Fo.sub.s] chart. Thus from the G values of the VS solution with the grout as the homogeneous medium, the underlying [theta] value is converted to the dimensionless G([k.sub.s][theta]/q') value, and the underlying time value is converted to the new Fp([[alpha].sub.s][tau]/[r*.sup.2]). This translation of the original plot essentially is rotation of the original G curve in the ratio of [k.sub.s]/[k.sub.g] followed by a parallel shift of this in the x-axis direction by Log([[alpha].sub.s]/[[alpha].sub.g]). This creates the "early stage" G value line. It is observed that the actual FEM results match with this early line fairly closely until about an Fo of 10 in this specific situation. The "late stage" line is obtained by adding a steady state differential borehole resistance ([DELTA][R.sub.b]) to the G value line for the step heat response in the homogeneous medium using the appropriate parameters in Equation 11. The [DELTA][R.sub.b] may be computed on the basis of multipole methodology (Hellstorm 1991). The example calculations are based on the approximate formula of Hellstrom (1991).

[FIGURE 10 OMITTED]

The differential borehole resistance ([DELTA][R.sub.b]) is the difference between the steady state thermal resistance of the borehole with the grout ([R.sub.b]) and the same with the soil replacing the grout, i.e., for the undistureed ground. This is positive when [k.sub.g] < [k.sub.s] and negative when [k.sub.g]>[k.sub.s] * [DELTA][R.sub.b] is added to the (G-curve response of the homogeneous medium essentially to model the process of replacing the medium in the borehole by a material with different properties (grout) in the drilling process. Within the range of values examined (detailed in Table 1), the computed [DELTA][R.sub.b] values compare reasonably well ([+ or -]6%) with the observed [DELTA][R.sub.b] values from the FEM results at around Fo = 100. Since this is only a differential value added to the homogeneous medium response, the effect of this error is reduced. Further, this error is systemic in nature; hence, an appropriate correction term may be added in future after a thorough parametric study.

The early line and late line intersect at a point corresponding to a transition time, as proposed by Sutton et al. (2002). When the grout and the medium have identical thermal properties (k,[alpha]), the early line and the late lines overlap, as the early line undergoes no rotation or shift. In the transition zone, to the left of the transition point, the G plot exits from the early line and merges with the late line. The late line corresponds to a period when the thermal front has reached the borehole boundary, and the thermal process is governed by the properties of the medium. The transition time depends on the borehole to tube diameter ratio and on other parameters. A larger ratio would clearly lead to a later transition. It is expected that the width of the transition band would depend on the ratio of the heat capacities of the grout and the medium, a higher grout heat capacity implying a wider band.

In Figure 11, a similar G plot is shown for identical conditions as for Figure 10, except that the grout thermal conductivity is doubled. It is seen that the early line and the late line intersects nearly at the same Fourier number of about 17-19, but in an absolute time scale, they represent widely differing values. While for the grout with the higher [k.sub.g] value, the transition time is at about 3 hours, the same for the other grout is about 14 hours. It is ascertained from Figures 10 and 11 that the early lines and late lines form the envelope for the G value in the transition zone, and a distinct transition time exists. The factors that determine the relative positions of the takeoff point and the merging point vis-a-vis the transition point on the early and the late lines, however, are still unclear and would require further parametric study. Similarly, the shape of the G plot and its analytic representation in the transition zone may be established in the future through further study.

[FIGURE 11 OMITTED]

CONCLUSION

In previous studies, the short time values of the water temperature in the GHX loop after applying a step heat pulse have been obtained through numerical methods. As these are complex procedures, their usefulness is limited for either routine design purposes or for system simulation. Further, all the currently used analytical solutions are based on steady flux assumption and ignore the aggregate thermal capacity of the fluid. An analytical solution developed for the measurement of thermal conductivity with a cylindrical probe has been applied successfully to derive the short time temperature response of the fluid with the U-tube geometry in the homogeneous media. The solution also incorporates the aggregate thermal capacity of the fluid and the effect of the separation of the two legs of the U-tube. A method for obtaining the water temperature in composite media is also presented. This method is capable of capturing the early-stage response and the late-stage response quite successfully and also points to a transition point. Since the expression of the analytical solution is in the form of a definite integral, it may be evaluated numerically facilitating easier integration in the design softwares or energy simulation packages. The computed design lengths are expected to be lower than corresponding values obtained using the classical solutions as this solution takes into account thermal capacity of the circulating fluid. Further research is required for accurately predicting the temperature behavior in the transition time zone and undertaking a comprehensive sensitivity analysis of the design lengths of the GHX for different type of buildings with differing load profiles.

ACKNOWLEDGMENT

The authors gratefully acknowledge the funding support received from the North Dakota Department of Commerce, Division of Community Services, State Energy Office for the studies. The authors are grateful to three anonymous referees for their valuable comments and suggestions.

NOMENCLATURE

t = temperature, [degrees]C ([degrees]F)

c = specific heat J/kg * K (BTU/lb * [degrees]R)

D = diameter, m (ft)

G = dimensionless temperature ((k[theta])/q) or dimensionless thermal resistance

h = convective heat transfer coefficient/film coefficient W/[m.sup.2] * K (Btu/h * [ft.sup.2][degrees]F)

J = Bessel function of the first kind

k = thermal conductivity W/m * K (Btu/h * ft * [degrees]F)

L = characteristic length m (ft)

LpR = loop ratio (ratio of total fluid mass in the system to the fluid mass in the boreholes)

q' = rate of heat generated (building load) per unit length of borehole W/m (Btu/h.ft)

r = radius m (ft)

R = thermal resistance [m.sup.2]*K/W ([ft.sup.2]*h*[degrees]F/BTU)

S = thermal capacity per unit length of a cylindrical body J/m*[degrees]C (BTU/ft*[degrees]F)

SpR = tube spacing ratio (ratio of distance between the two tube walls to the tube outer radius)

u = integration variable (dummy variable)

Y = Bessel function of the second kind

Dimensionless Numbers

Bi = Biot number (hL/k)

Fo = Fourier number ([alpha]t/[r*.sup.2]

Greek Symbols

[alpha] = thermal diffusivity [ft.sup.2]/h ([m.sup.2]/s)

[gamma] = Euler's constant (0.5772 ...)

[theta] = (t-[t.sub.[infinity]]): temperature excess (above or below the undistributed ground temperature)

[rho] = density kg/[m.sup.3] (lb/[ft.sup.3])

[tau] = time, s

Subscripts/Superscripts

* = core

b = borehole

bb = borehole boundary

eq = equivalent

f = fluid

g = grout (filling in the borehole)

i = inside

o = outside

p = pipe

s = soil/medium (ground)

sf = steady flux

REFERENCES

Blackwell, J.H. 1954. A Transient-Flow method for determination of thermal constants of insulating materials in bulk. Journal of Applied Physics. 25(2): 137-44.

Bose, J.H. 1984. Closed-loop/ Ground source heat pump design manual. Engineering Technology extension, Oklahoma State University Stillwater, OK.

Carslaw, H. S., & Jaeger, J. C. (1959). Conduction of heat in solids. Oxford: Clarendon Press.

Claesson, J., Dunard, A.1983. Heat extraction from the ground by horizontal pipes- a mathematical analysis. Swedish Council of Building Research, document D1, Stockholm, Sweden.

Deerman, J. D. and Kavanaugh, S. P. 1991. Simulation of vertical u-tube ground coupled heat pump systems using the cylindrical heat source solution. ASHRAE Transactions 97(1): 287-295.

Gu, Yian, and O'Neal, D.L.1995. An analytical solution to transient heat conduction in a composite region with a cylindrical heat source. Journal of Solar Energy Engineering, Transactions of the ASME 117 (3):242-248.

Gu, Y., and D. O'Neal. 1998. Development of an Equivalent Diameter Expression for Vertical U-Tubes Used in Ground-Coupled Heat Pumps. ASHRAE Transactions 104(2):347-355.

Gehlin, S.E.A., Hellstrom, G. 2003 Comparison of four models for thermal response test evaluation. ASHRAE Transactions 109(2):131-142.

Hellstrom, G. (1991). Ground Heat Storage. Thermal Analyses of Duct Storage Systems I: Theory. University of Lund, Department of Mathematical Physics. Lund, Sweden.

Kavanaugh, S.P. (1984). Simulation and experimental verification of vertical ground-coupled heat pump systems, Ph.D Dissertation, Oklahoma State University, Stillwater, OK.

Kavanaugh, S. P., K. Rafferty. 1997. Ground Source Heat Pumps: Design of Geothermal Systems for Commercial and Institutional Buildings. Atlanta: American Society of Heating, Refrigerating and Air-Conditioning Engineers, Inc.

Lamarche, L., Beauchamp, B. 2007. New solutions for the short-time analysis of geothermal vertical boreholes. International Journal of Heat and Mass Transfer 50 (7-8):1408-1419.

Myers, G.E.1998.Analytical Methods in Conduction Heat Transfer. AMCHT publications, Madison, WI.

Mei, V. C., Baxter, V. D., 1986. Performance of a Ground Coupled Heat Pump with Multiple Dissimilar U-Tube Coils in Series. ASHRAE Transactions, 92(2A):30-42.

Muraya, N. K. (1994). Numerical modeling of the transient thermal interference of vertical U-tube heat exchangers. Ph.D. Dissertation. Texas A&M University, College Station, TX.

Rottmayer, S. P., W. A. Beckman, and J. W. Mitchell. (1997). Simulation of a single vertical U-tube ground heat exchanger in an infinite medium. ASHRAE Transactions. 103(2):651-659.

Shonder J.A., J. V. Beck. 1999. Determining Effective Soil Formation Thermal Properties from Field Data Using Parameter Estimation Technique. ASHRAE Transactions 105(1).

Sutton, M.G., Couvillion, R.J., Nutter, D.W., Davis, R.K. 2002. An algorithm for approximating the performance of vertical bore heat exchangers installed in a stratified geological regime. ASHRAE Transactions 108(2): 177-184.

Xu, X., Spitler, J.D., Modeling of Vertical Ground Loop Heat Exchangers With Variable Convective Resistance And Thermal Mass of The Fluid The Tenth International Conference on Thermal Energy Storage (ECOSTOCK), Richard Stockton College, Pomona, NJ, May 31 - June 2, 2006.

Yavuzturk, C., and J. D. Spitler. (1999). A Short Time Step Response Factor Model for Vertical Ground Loop Heat Exchangers. ASHRAE Transactions. 105(2):475-485.

Young, R. (2004). Development, Verification, and Design Analysis of the Borehole Fluid Thermal Mass Model for Approximating Short Term Borehole Thermal Response. Master's Thesis. Oklahoma. State Univ. Stillwater, OK.

Gopal Bandyopadhyay, PhD

Associate Member ASHRAE

Manohar Kulkarni, PhD, PE

Member ASHRAE

Michael Mann, PhD

Gopal Bandyopadhyay is a postdoctoral research associate and Michael Mann is professor and department chair of the Department of Chemical Engineering and Manohar Kulkarni is professor and department chair of the Department of Mechanical Engineering at the University of North Dakota, Grand Forks, ND.

Printer friendly Cite/link Email Feedback | |

Author: | Bandyopadhyay, Gopal; Kulkarni, Manohar; Mann, Michael |
---|---|

Publication: | ASHRAE Transactions |

Date: | Jul 1, 2008 |

Words: | 8867 |

Previous Article: | Emulating nature: evaporative cooling systems. |

Next Article: | A numerical estimate of flexible short-tube flow and deformation with R-134a and R-410a. |