

small (250x250 max)
medium (500x500 max)
Large
Extra Large
large ( > 500x500)
Full Resolution


DEVELOPMENT, VERIFICATION, AND DESIGN ANALYSIS OF THE BOREHOLE FLUID THERMAL MASS MODEL FOR APPROXIMATING SHORT TERM BOREHOLE THERMAL RESPONSE by THOMAS RAY YOUNG Bachelor of Science Oklahoma State University Stillwater, Oklahoma 2001 Submitted to the Faculty of the Graduate College of the Oklahoma State University in partial fulfillment of the requirements for the degree of MASTERS OF SCIENCE Oklahoma State University December, 2004 ii DEVELOPMENT, VERIFICATION, AND DESIGN ANALYSIS OF THE BOREHOLE FLUID THERMAL MASS MODEL FOR APPROXIMATING SHORT TERM BOREHOLE THERMAL RESPONSE Thesis Approved: Dr. Spitler Thesis Advisor Dr. Delahoussaye Dr. Fisher Dr. Emsile Dean of Graduate College iii ACKNOWLEDGEMENTS In order to give credit where it is due, my conscience requires that, at the very least, I mention God who sent his son Jesus Christ to die in my place, for my sins, and gives eternal life to everyone who believes in Jesus. Without the grace, mercy, strength and wisdom God has provided me in various ways, this thesis would never have been finished. He deserves ALL the honor and glory. Also I am indebted to Dr. Spitler for his willingness to take me on as one of his students. His expertise and experience were invaluable and his willingness to continue working with me after I took a job and left Stillwater is very much appreciated I would like to thank Dr. Delahoussaye for his mentoring early on in my research. Dr. Delahoussaye was essential in helping me gain the programming skills and practical understanding that I needed to succeed. To Aditya, Hayder, and Liu I appreciate your friendships, learned from your expertise, and enjoyed working with you greatly. I would especially like to thank Xiaowei for the hours he spent running simulation and examining my work. I would also like to thank my wife Rachel who provided moral support through a listening ear and plenty of cookies and brownies. I would like to thank my parents for instilling in me moral values and a good work ethic. Also, thank you for your prayers and financial support. iv I should also mention my Sunday school class for all there prayers and moral support. Hardly a Sunday went by without someone asking me how the thesis was going. v TABLE OF CONTENTS 1 INTRODUCTION ....................................................................................................1 1.1 BACKGROUND....................................................................................................................... 5 1.2 LITERATURE REVIEW......................................................................................................... 5 1.2.1 Steady State Modeling of Boreholes ....................................................................................... 7 1.2.1.1. Line Source Model..................................................................................................................... 11 1.2.1.2. GuO’Neal Equivalent Diameter Model .................................................................................... 15 1.2.1.3. Paul Model................................................................................................................................ 17 1.2.1.4. Cylinder Source Model .............................................................................................................. 19 1.2.1.5. Multipole................................................................................................................................... 21 1.2.2 Transient Modeling of Borehole Heat Exchangers............................................................... 31 1.2.2.1. Buried Electrical Cable Model................................................................................................... 32 1.2.2.2. Gfunction Model: Long Time Step........................................................................................... 35 1.2.2.3. Gfunction Model: Short Time Step........................................................................................... 41 1.2.3 GLHEPRO Version 3 Design Tool ....................................................................................... 47 1.3 OBJECTIVES......................................................................................................................... 48 2 COMPARISON OF BOREHOLE RESISTANCE CALCULATION METHODS...................................................................................................................... 50 2.1 BOREHOLE RESISTANCE TRANSIENT AND STEADY STATE ..................................................... 51 2.2 BOREHOLE RESISTANCE CALCULATION FROM ANALYTICAL AND EMPIRICAL METHODS... 53 2.3 BOREHOLE RESISTANCE CALCULATION USING NUMERICAL METHODS ............................... 59 2.4 NUMERICAL METHODS: COMPARISON BETWEEN GEMS2D AND THE PIESECTOR APPROXIMATION FOR CALCULATING STEADY STATE RESISTANCE..................................................... 59 2.5 COMPARISON OF METHODS FOR CALCULATING STEADY STATE BOREHOLE RESISTANCE.. 62 2.6 BOREHOLE RESISTANCE AND MERGING OF THE SHORT AND LONG TIME STEP GFUNCTION 68 2.7 CONCLUSION............................................................................................................................ 71 3 SHORT TIME STEP GFUNCTION CREATION AND THE BOREHOLE FLUID THERMAL MASS MODEL (BFTM)............................................................. 73 3.1 BOREHOLE FLUID THERMAL MASS MODEL ........................................................................... 74 3.2 GROUT ALLOCATION FACTOR USED TO IMPROVE ACCURACY ............................................. 77 3.3 FLUID MULTIPLICATION FACTOR IN THE BFTM MODEL....................................................... 78 3.4 IMPLEMENTATION OF THE BFTM MODEL.............................................................................. 79 3.4.1 Bessel Function Evaluation .................................................................................................. 80 3.4.2 BFTM Model  Solving the Integral...................................................................................... 81 3.4.3 Incorporating the Fluid Thermal Mass Model in a Design Program................................... 83 3.5 IMPROVING THE BFTM MODEL FOR SMALL TIMES USING LOGARITHMIC EXTRAPOLATION 86 3.5.1 Implementing Logarithmic Extrapolation............................................................................. 87 4 NUMERICAL VALIDATION OF THE BOREHOLE FLUID THERMAL MASS MODEL USING GEMS2D................................................................................ 90 vi 4.1 GEMS2D SIMULATIONS .......................................................................................................... 93 4.2 FINDING THE GROUT ALLOCATION FACTOR .......................................................................... 95 4.3 BOREHOLE DIAMETER VALIDATION WITH LINE SOURCE COMPARISON ............................ 100 4.4 SHANK SPACING VALIDATION ............................................................................................... 105 4.5 GROUT CONDUCTIVITY VALIDATION.................................................................................... 110 4.6 SOIL CONDUCTIVITY VALIDATION ........................................................................................ 112 4.7 GROUT VOLUMETRIC HEAT CAPACITY VALIDATION .......................................................... 116 4.8 BFTM MODEL FLUID FACTOR VALIDATION WITH GEMS2D............................................. 120 4.9 IMPLEMENTATION AND VALIDATION OF THE BFTME MODEL.......................................... 124 4.10 CONCLUSION OF BFTM MODEL VALIDATION...................................................................... 129 5 THE EFFECT OF THE BFTM MODEL ON GLHE DESIGN....................... 130 5.1 TEST BUILDINGS..................................................................................................................... 130 5.1.1 Church ............................................................................................................................... 130 5.1.2 Small Office Building.......................................................................................................... 132 5.1.3 Annual Loading .................................................................................................................. 133 5.2 GLHE DESIGN PROCEDURES................................................................................................. 141 5.3 SIMULATION RESULTS............................................................................................................ 147 5.3.1 Fluid Factor Results ........................................................................................................... 148 5.3.2 UTube Shank Spacing Results........................................................................................... 159 5.3.3 Borehole Diameter Results ................................................................................................. 164 5.1.4 Discussion of GLHEPRO Sizing Results ............................................................................ 168 6 HOURLY SIMULATION USING THE BFTM MODEL................................ 170 6.1 HVACSIM+ HOURLY SIMULATION...................................................................................... 171 6.2 LINE SOURCE AND BFTM MODEL COMPARISON USING A DETAILED HVACSIM+ MODEL 173 6.3 INFLUENCE OF THE FLUID MULTIPLICATION FACTOR ON SYSTEM DESIGN........................ 181 6.3.1 Fluid Factor Analysis with HVACSIM+ Simulation Tools................................................. 182 6.4 HVACSIM+ AND GLHEPRO COMPARISON........................................................................ 185 6.5 CONCLUSION.......................................................................................................................... 193 7 CONCLUSIONS AND RECOMMENDATIONS.............................................. 194 7.1 CONCLUSIONS........................................................................................................................ 194 7.2 RECOMMENDATIONS .............................................................................................................. 196 vii LIST OF TABLES Table Page Table 11 Paul Curve Fit Parameters used to Calculate the Steady State Grout Resistance .................................................................................................................................. 18 Table 12 Variable Input List for the Multipole Method.................................................. 23 Table 21 Borehole Properties (Base Case) ...................................................................... 60 Table 22 Borehole Resistance Comparison between GEMS2D and the pie sector approximation ........................................................................................................... 61 Table 23 Base Line Borehole Properties ......................................................................... 63 Table 24 Steady State Borehole Resistance Comparison................................................ 65 Table 25 Percent Error of Borehole Resistance............................................................... 66 Table 31 Borehole Properties Table ................................................................................ 75 Table 41 Borehole Properties for GEMS2D to BFTM Model Comparisons. ................. 91 Table 42 Borehole Diameter vs Time for Slope Matching. ............................................ 98 Table 43 GAF Dependent on Borehole Diameter, Shank Spacing and Fluid Factor...... 99 Table 44 Borehole Resistances and GAF for Diameter Validation Tests ..................... 101 Table 45 Borehole Resistances for Shank Spacing Validation Tests............................ 106 Table 46 Borehole Resistances for Shank Spacing Validation Tests............................ 110 Table 47 Borehole Resistances for Soil Conductivity Validation Tests ....................... 113 Table 48 NonDimensional Borehole Diameter vs Time for Slope Matching.............. 126 Table 49 GAF Dependent on NonDimensional Borehole Diameter, NonDimensional Shank Spacing and Fluid Factor ............................................................................. 126 Table 51 Church Building Description.......................................................................... 131 Table 52 Church Building Load Table for Different Locations.................................... 134 Table 53 Building Load Table for the Small Office Building....................................... 138 Table 54 GLHE Properties for the Church and Small Office Building......................... 141 Table 55 Undisturbed Ground Temperature Table for Various Cities.......................... 142 Table 56 Percent Change in GLHEPRO Sizing Depth with Respect to Varying Fluid Factor for PeakLoadDominant and NonPeakLoadDominant Buildings with Standard Grout ........................................................................................................ 154 Table 57 Percent Change in GLHEPRO Sizing Depth with Respect to Varying Fluid Factor for PeakLoadDominant and NonPeakLoadDominant Buildings with Thermally Enhanced Grout..................................................................................... 155 Table 58 Required Depth (ft) for Thermally Enhanced Grout, Calculated with GLHEPRO.............................................................................................................. 156 Table 59 Required Depth (ft) for Standard Grout, Calculated with GLHEPRO........... 157 viii Table 510 Percent Change in GLHEPRO Sizing Depth for Changes in Borehole Shank Spacing.................................................................................................................... 163 Table 511 Percent Change in GLHEPRO Sizing Depth for Changes in Borehole Diameter.................................................................................................................. 167 Table 61 Coefficients for the VS200 Climate Master Heat Pump ................................ 172 ix LIST OF FIGURES Figure Page Figure 11 Borehole system ................................................................................................ 2 Figure 12 Borehole Cross Section..................................................................................... 3 Figure 13 Crosssection of the Borehole and the Corresponding Thermal ΔCircuit (Hellström, 1991 p.78)................................................................................................ 7 Figure 14 Cross Section of a Borehole with Symmetry Line and the Corresponding Thermal Circuit........................................................................................................... 8 Figure 15 Actual Geometry vs Equivalent Diameter Approximation............................. 16 Figure 16 Types of shank spacing used in the Paul borehole resistance approximation. 18 Figure 17 Example of a 2D System for the Multipole Method. ...................................... 22 Figure 18 Source and Sink Location for a Single Pipe.................................................... 26 Figure 19 Steady State Temperature Field for a Borehole Heat Exchanger ................... 29 Figure 110 Temperature Change Around the Borehole Circumference.......................... 30 Figure 111 Diagram of a buried electrical cable and circuit ........................................... 33 Figure 112 Short and long time step gfunction without borehole resistance ................. 37 Figure 113 Short and long time step gfunctions with borehole resistance..................... 37 Figure 114 Twodimensional radialaxial mesh for a heat extraction borehole in the ground (Eskilson, 1987)............................................................................................ 39 Figure 115 Long Time Step gfunction for a 64 Borehole System in an 8x8 Configuration with Varying Borehole Spacing ........................................................ 41 Figure 116 Grid for a cross section of a borehole ........................................................... 43 Figure 117 Grid for PieSector Approximation (Yavuzturk, 1999) ................................ 45 Figure 21 Transient Borehole Resistance Profile vs Time .............................................. 52 Figure 22 Percent Difference in Transient Borehole Resistance with respect to Steady State borehole resistance........................................................................................... 52 Figure 23 Cylinder Source Diagram for Calculating the Utube Outside Wall Temperature for use in the Steady State Borehole Resistance Calculation .............. 55 Figure 24 Line Source STS Gfunction Compared to LTS Gfunction Using Different Borehole Resistance Calculation Methods for a Single Borehole System ............... 70 Figure 31 Average Fluid Temperature using the line source and GEMS2D model with fluid mass for a heat rejection pulse ......................................................................... 73 Figure 32 Integrated Function in the BEC Model ........................................................... 83 Figure 33 STS GFunction Translation ........................................................................... 85 Figure 34 BFTME, BFTM, and GEMS2D Fluid Temperature and Gfunction ............ 89 Figure 41 Borehole Geometries Simulated with GEMS2D ............................................ 92 Figure 42 Blocks for a Borehole System without Interior Cells for GEMS2D............... 93 x Figure 43 GEMS2D Grid of a Borehole with Soil .......................................................... 94 Figure 44 GEMS2D Grid of Utube and Fluid with 8 Annular Regions ........................ 95 Figure 45 Fluid Temperature vs Time of GEMS2D and BFTM with Varying GAF...... 96 Figure 46 Gfunction Slope vs GAF for the BFTM at 5 Hours ...................................... 97 Figure 47 Fluid Temperature and Gfunction for 15.24 cm (6 in) Diameter Borehole Using a Slope Matching Time of 3 Hours for (a) and (b) and 5 hours for (c) and (d) .................................................................................................................................. 98 Figure 48 Validation of the BFTM Model Using a GEMS2D Simulation with a 7.62 cm (3 in) Borehole ........................................................................................................ 102 Figure 49 Validation of the BFTM Model Using a GEMS2D Simulation with a 11.4 cm (4.5 in) Borehole ..................................................................................................... 103 Figure 410 Validation of the BFTM Model Using a GEMS2D Simulation with 15.2 cm (6 in) Borehole with 3.16 cm (1.24 in) Shank Spacing .......................................... 104 Figure 411 Validation of the BFTM Model Using a GEMS2D Simulation with a 19.1 cm (7.5 in) Borehole with a 4.12 cm (1.62 in) Shank Spacing..................................... 105 Figure 412 Validation of the BFTM Model Using a GEMS2D Simulation with a 0.316 cm (0.125 in) Shank Spacing ...................................................................................107 Figure 413 Validation of the BFTM Model Using a GEMS2D Simulation with a 2.25 cm (0.89 in) Shank Spacing .................................................................................... 108 Figure 414 Validation of the BFTM Model Using a GEMS2D Simulation with a 1 4.13 cm (5/8 in) Shank Spacing...................................................................................... 109 Figure 415 Validation of the BFTM Model Using a GEMS2D Simulation with Grout Conductivity of 0.25 W/(m·K) (0.144 Btu/(h·ft·°F)) .............................................. 111 Figure 416 Validation of the BFTM Model Using a GEMS2D Simulation with grout conductivity of 1.5 W/m·K (0.867 Btu/(h·ft·°F)).................................................... 112 Figure 417 Validation of the BFTM Model Using a GEMS2D Simulation With Soil Conductivity of 0.5 W/m·K (0.289 Btu/(h·ft·°F))................................................... 114 Figure 418 Validation of the BFTM Model Using a GEMS2D Simulation with a Soil Conductivity of 1.5 W/m·K (0.867 Btu/(h·ft·°F))................................................... 115 Figure 419 Validation of the BFTM Model Using a GEMS2D Simulation with a Soil Conductivity of 8 W/m·K (4.62 Btu/(h·ft·°F))........................................................ 116 Figure 420 Validation of the BFTM Model Using a GEMS2D Simulation with a Grout Volumetric Heat Capacity of 2 MJ/m3·k (29.8 Btu/ft3·F) ....................................... 118 Figure 421 Validation of the BFTM Model Using a GEMS2D Simulation with a Grout Volumetric Heat Capacity of 8 MJ/m3K (119 Btu/ft3·F)....................................... 119 Figure 422 Validation of the BFTM Model Using a GEMS2D Simulation with a 11.4 cm (4.5 in) BH Diameter 3 cm (1.18 in) Shank Spacing and 2 Times the Fluid.......... 121 Figure 423 Validation of the BFTM Model Using a GEMS2D Simulation with a 11.4 cm (4.5 in) BH Diameter 3 cm (1.18 in) Shank Spacing and 4 Times the Fluid.......... 122 Figure 424 Validation of the BFTM Model Using a GEMS2D Simulation with a 19.05 cm (7.5 in) BH Diameter 2.25 cm (0.886 in) Shank Spacing and 2 Times the Fluid ................................................................................................................................ 123 Figure 425 Validation of the BFTM Model Using a GEMS2D Simulation with a 19.05 cm (7.5 in) BH Diameter 2.25 cm (0.886 in) Shank Spacing and 4 Times the Fluid ................................................................................................................................ 124 xi Figure 426 Temperature Profile for the BFTME and GEMS2D Models with a 17.8 cm (7 in) Borehole Diameter and 4 cm (1.57 in) Shank Spacing Using an Interpolated GAF Value.............................................................................................................. 128 Figure 51 Monthly Church Heating Loads.................................................................... 135 Figure 52 Monthly Church Cooling Loads.................................................................... 136 Figure 53 Monthly Church Peak Heating Loads ........................................................... 137 Figure 54 Monthly Church Peak Cooling Loads........................................................... 137 Figure 55 Monthly Heating Loads for the Small Office Building ................................ 139 Figure 56 Monthly Cooling Loads for the Small Office Building ................................ 139 Figure 57 Monthly Peak Heating Loads for the Small Office Building........................ 140 Figure 58 Monthly Peak Cooling Loads for the Small Office Building ....................... 140 Figure 59 Raw Church Loads Single 2 Hour Peak Heat Load...................................... 144 Figure 510 Raw Church Loads for Birmingham AL..................................................... 144 Figure 511 One Peak of Hourly Loads for Tulsa Small Office Building...................... 145 Figure 512 One Work Week of Hourly Loads for Tulsa Small Office Building .......... 146 Figure 513 Short Time Step GFunction Comparison between the Line Source and the BFTM model with 0.1, 1, 2, 3 and 4 x Fluid Factor ............................................... 147 Figure 514 GLHEPRO Sized Depth vs Fluid Factor for a Church Building with a Borehole Diameter of 4.5 in (11.4 cm), Enhanced Grout, and Shank Spacing of 0.125 in (0.318 cm)................................................................................................. 149 Figure 515 GLHEPRO Sized Depth vs Fluid Factor for Church Building with a Borehole Diameter of 4.5 in (11.4 cm), Enhanced Grout, and Shank Spacing of 1.87 in (4.75 cm)............................................................................................................. 149 Figure 516 GLHEPRO Sized Depth vs Fluid Factor for a Church Building with a Borehole Diameter of 4.5 in (11.4 cm), Standard Bentonite Grout, and Shank Spacing of 0.125 in (0.318 cm)............................................................................... 150 Figure 517 GLHEPRO Sized Depth vs Fluid Factor for Church Building with a Borehole Diameter of 4.5 in (11.4 cm), Standard Bentonite Grout, and Shank Spacing of 1.87 in (4.75 cm)................................................................................... 151 Figure 518 GLHEPRO Sized Depth vs Fluid Factor for Small Office Building with a Borehole Diameter of 4.5 in (11.4 cm), Standard Bentonite Grout, and Shank Spacing of 0.125 in (0.318 cm)............................................................................... 152 Figure 519 GLHEPRO Sized Depth vs Fluid Factor for Small Office Building with a Borehole Diameter of 4.5 in (11.4 cm), Thermally Enhanced Grout, and Shank Spacing of 1.87 in (4.75 cm)................................................................................... 152 Figure 520 Borehole Resistance vs Shank Spacing for Church Building Using Standard Bentonite Grout and a Diameter of 4.5 in (11.4 cm) .............................................. 159 Figure 521 Sized Borehole Depth vs Shank Spacing for Church Building Using Standard Bentonite Grout, Diameter of 4.5 in (11.4 cm), Fluid Factor of 1.......................... 160 Figure 522 Sized Borehole Depth vs Shank Spacing for Church Building Using Standard Bentonite Grout, Diameter of 4.5 in (11.4 cm), Fluid Factor of 2.......................... 161 Figure 523 Sized Borehole Depth vs Shank Spacing for Office Building Using Standard Bentonite Grout, Diameter of 4.5 in (11.4 cm), Fluid Factor of 1.......................... 161 Figure 524 Sized Borehole Depth vs Shank Spacing for Office Building Using Standard Bentonite Grout, Diameter of 4.5 in (11.4 cm), Fluid Factor of 2.......................... 162 xii Figure 525 Borehole Resistance vs Diameter Using Thermally Enhanced Grout, Shank Spacing of 0.125 in (0.318 cm)............................................................................... 164 Figure 526 Sized Borehole Depth vs Diameter for Church Building Using Thermally Enhanced Grout, Shank Spacing of 0.125 in (0.318 cm), and Fluid Factor of 1.... 165 Figure 527 Sized Borehole Depth vs Diameter for Small Office Building Using Standard Bentonite Grout , Shank Spacing of 0.125 in (0.318 cm), and Fluid Factor of 1 ... 165 Figure 61 Three Component Model of a GLHE System for Hourly Loads .................. 171 Figure 62 Gfunction’s for Various Fluid Factors......................................................... 173 Figure 63 Detailed HVACSIM+ Model with Tulsa Loads for PeakLoadDominant Times....................................................................................................................... 175 Figure 64 Detailed HVACSIM+ Model with Tulsa Loads for PeakLoadDominant Times....................................................................................................................... 175 Figure 65 Detailed HVACSIM+ with Tulsa Loads....................................................... 176 Figure 66 Detailed HVACSIM+ with Tulsa Loads for a Short Duration Heat Pulse on a Long Duration Heat Pulse....................................................................................... 177 Figure 67 Detailed HVACSIM+ with Tulsa Loads for Long Duration Heat Pulses..... 178 Figure 68 Difference in Temperature between the LS and BFTM Models................... 179 Figure 69 Heat Pump Power Curve for the LS and BFTM Models.............................. 180 Figure 610 Heat Pump Power Curve for the LS and BFTM Models............................ 180 Figure 611 Detailed HVACSIM+ Model with Tulsa Loads ......................................... 181 Figure 612 GLHE Inlet Temperature for Different Fluid Factors................................. 183 Figure 613 GLHE Inlet Temperature for Different Fluid Factor .................................. 184 Figure 614 GLHE Inlet Temperature for Different Fluid Factors................................. 185 Figure 615 Typical Peak Loads for Small Office Building in Houston ........................ 186 Figure 616 Typical Peak Loads for Church Building in Nashville ............................... 187 Figure 617 Two Component GLHE Model................................................................... 187 Figure 618 Entering Temperature to the Heat Pump for a Church Building Located in Nashville, 1 x fluid, and 2 hour GLHEPRO Peak Duration ................................... 188 Figure 619 Entering Temperature to the Heat Pump for a Church Building Located in Nashville, 2 x fluid, and 2 hour GLHEPRO Peak Duration ................................... 189 Figure 620 Entering Temperature to the Heat Pump for a Church Building Located in Nashville, 4 x fluid, and 2 hour GLHEPRO Peak Duration ................................... 189 Figure 621 GLHEPRO and HVACSIM+ Maximum Yearly Entering Temperature to the Heat Pump for a Church Building Located in Nashville ........................................ 190 Figure 622 GLHEPRO and HVACSIM+ Minimum Yearly Entering Temperature to the Heat Pump for a Church Building Located in Nashville ........................................ 191 Figure 623 Entering Temperature to the Heat Pump for a Small Office Building Located in Houston, 1 x fluid, and 8 hour GLHEPRO Peak Duration................................. 192 Figure 624 Entering Temperature to the Heat Pump for a Small Office Building Located in Houston, 2 x fluid, and 8 hour GLHEPRO Peak Duration................................. 192 1 1 INTRODUCTION “As long as the earth endures, seedtime and harvest, cold and heat, summer and winter, day and night will never cease.” (Genesis 8:22, NIV) The basic needs of man which include keeping warm in winter and cool in summer have remained constant throughout history. The technology used to meet the needs has changed. People and animals have historically used caves and manmade holes as shelter from the elements. In this way humans have been extracting heat from the earth to keep warm in winter and using the earth to keep cool in summer for centuries. Modern man uses more refined methods for extracting and rejecting heat from the ground such as ground source heat pump (GSHP) systems. The term “ground source heat pump (GSHP)” refers to heat pump systems that use either the ground or a water reservoir as a heat source or sink. GSHP systems are either openloop or closedloop. Openloop GSHP systems use a pump to circulate groundwater through the heat pump heat exchanger. A closedloop GSHP system uses a water pump to circulate fluid through pipes buried horizontally or inserted into boreholes in the ground. The buried closed loop version of the GSHP is commonly referred to as a ground loop heat exchanger (GLHE). The physical properties of boreholes are very important to the study of GLHE systems. Boreholes typically range between 46 to 122 meters (150 to 400 ft) deep and are typically around 10 to 15 cm (4 to 6 inches) in diameter. A borehole system can be 2 composed of anywhere from 1 to over 100 boreholes. Each borehole in a multiborehole system is typically placed at least 4.5 m (15 ft) from all other boreholes. Figure 11 shows a vertical cross section of three boreholes. Each borehole is connected to the other boreholes with pipes that are typically buried 12 meters (36 ft) under the top surface. Figure 11 Borehole system After the Utube is inserted, the borehole will usually be backfilled with grout. The grout is used to prevent contamination of aquifers. Figure 11 shows 3 ideal boreholes. The grout is often a bentonite clay mixture, with the possibility of having thermallyenhanced additives. The grout usually has a thermal conductivity significantly lower than the surrounding ground. The circulating fluid is water or a waterantifreeze mixture. Each borehole is shown in this picture to be parallel with the other boreholes and perpendicular to the surface of the ground. In reality, drilling rigs do not drill perfectly straight, causing the path of a borehole to deviate, especially in deep boreholes. The Utube as shown in Figure 11 has equal spacing between the two legs of the Utube throughout the borehole. In real systems, however, the Utube leg spacing does 3 not necessarily remain constant throughout the length of the borehole. Spacers are sometimes employed to force the tubes towards the borehole wall. Figure 12 shows a two dimensional horizontal cross section of a single borehole. The Utube leg spacing is called the shank spacing and is defined as the shortest distance between the outer pipe walls of each leg of the Utube. UTUBE SHANK SPACING GROUT BOREHOLE GROUND Figure 12 Borehole Cross Section As previously mentioned, the size of a borehole heat exchanger system can range from one borehole to over a hundred. For small buildings one borehole may suffice but for large commercial buildings over 100 boreholes are sometimes required. This can make the initial investment quite costly. The main advantage of a GSHP system over an airsource heat pump system is that it rejects heat to the ground in the summer, when the ground is cooler than the air, and extracts heat from the ground during the winter, when the ground is warmer than the air. A GSHP system will very seldom reject the same amount of heat as it extracts on an annual basis. In cold climates, for envelopedominated buildings, the GSHP system will extract much more heat from the ground than it rejects to the ground. In this case, the ground surrounding the boreholes gradually declines in temperature. Over time, the reduction in the ground temperature around the boreholes will decrease the performance 4 of the heat pump in heating mode. In cold climates, the fluid circulating in the boreholes might drop below freezing, requiring the addition of antifreeze in the system. Similarly in warm climates, since more heat is rejected to the ground than extracted from the ground, the ground temperature will rise. This will impair the performance of the heat pump in cooling mode. The actual annual imbalance depends not only on climate but also on the building internal heat gains and building design. The thermal loads over a number of years must be accounted for when designing a GSHP system. This is necessary to determine the impact of any annual heat imbalance. If a borehole heat exchanger (BHE) is overdesigned, the initial construction costs may be excessive. If the system is underdesigned, the BHE may not meet the long term heating or cooling needs of the user. Research has been conducted for the purpose of applying GSHP technology to other areas besides buildings. Chiasson and Spitler (2001 and 2000) at Oklahoma State University have conducted research applying GSHP technology to highway bridges. The system uses pipes embedded in the road pavement to circulate fluid from a GSHP to eliminate ice or snow formation. The potential benefits of this new application include safer driving conditions and longer lasting bridges and roads due to reduced corrosion. Engineers who are attempting to design a GSHP system for a specific application can use programs such as GLHEPRO (Spitler 2000) to describe a potential system composed of a specific ground loop heat exchanger and heat pump and then simulate the systems response to monthly and peak, heating and cooling loads. Using programs such as GLHEPRO, engineers can also optimize the depth of a specific borehole heat exchanger configuration. 5 1.1 BACKGROUND Regardless of the GSHP application, the thermal response of the GLHE plays an important part in the design and simulation of GSHP systems. Since thermal loading on typical GLHE systems is of long duration, design methodologies have focused, in some detail, on long time step responses to monthly loads. (Eskilson, 1987) However, short time responses have typically been modeled only crudely using analytical models such as the cylinder source. These short time responses can be very important in determining the effect of daily peak loads. Daily peak loads occur in all applications but may be dominant in applications such as church buildings, concert halls, and the Smart Bridge application. To model the shorttime response, it is important to accurately represent such details as the borehole radius, Utube diameter and shank spacing, as well as the thermal properties and mass of the circulating fluid, Utube and grout. This thesis presents a new methodology for modeling the short time GLHE thermal response. This is particularly important for systems with peakloaddominant loading conditions. 1.2 LITERATURE REVIEW The literature review describes different methods that have been developed to model borehole heat exchangers. The methods are divided into two categories: steady state and transient. Quasisteady state conditions occur in twodimensional borehole cross sections, as shown in Figure 12 when the circulating fluid, Utubes and grout within a borehole do not change temperature (relative to each other) with time for a constant heat flux. If the internal borehole temperature differences are constant, the borehole resistance, defined as the resistance between the circulating fluid and the borehole, is also constant. Thus, 6 when a borehole’s internal temperature differences have stabilized for a constant heat flux, the borehole resistance can be modeled as a constant. Transient modeling of borehole heat exchangers might be broken into three different regions. The first region deals with transients that occur within the borehole before the borehole reaches steady state. For this transient region, the borehole may be modeled as having infinite length since surface and bottom end effects can be neglected. Two dimensional geometric and thermal properties of the borehole influence the temperature response in the first region. The second region occurs after the internal geometric and thermal properties of the borehole cease to influence the temperature response and before the surface and bottom end effects influence the temperature response. The third transient region occurs when three dimensional effects such as borehole to borehole interaction, surface and bottom end effects influence the temperature response. The borehole transient resistance or gfunction is broken into two zones called the short time step (STS) gfunction (Yavuzturk, 1999) and the long time step (LTS) gfunction. (Eskilson, 1987) The short and the long time step gfunctions relate to the three regions described above in that the short time step gfunction represents region one and two and the long time step gfunction represents region two and three. Thus it is important to note that the short and long time step gfunction can both represent region 2. This allows the two gfunctions to be integrated into one continuous gfunction curve, allowing the borehole transient resistance to be known for small times, such as 0.5 hours, to large times, such as 100 years. Short and long time step gfunctions are discussed in 7 detail in section 1.2.2.3 and 1.2.2.2 respectively. The next two sections, however, will discuss the current literature for steady state and transient borehole modeling. 1.2.1 Steady State Modeling of Boreholes This section discusses borehole resistance since it is an important part of transient analysis. The borehole resistance is the thermal resistance between the fluid and the borehole wall. Figure 13 shows a crosssection of a borehole and a corresponding thermal delta circuit. Figure 13 Crosssection of the Borehole and the Corresponding Thermal ΔCircuit (Hellström, 1991 p.78) f 1 T and f 2 T (°C or °F) represent the fluid temperature in each leg of the Utube and 1 q and 2 q ⎥⎦ ⎤ ⎢⎣ ⎡ hr ⋅ ft or Btu m W the heat flux (heat transfer rate per unit length of borehole) from the circulating fluid. b T represents the average temperature on the borehole wall. As shown in the delta circuit in Figure 13, the thermal resistance between f 1 T and b T is 1 R ⎥⎦ ⎤ ⎢⎣ ⎡ ⋅ ⋅° Btu or hr ft F W mK and the thermal resistance between f 2 T and b T is 2 R 8 ⎥⎦ ⎤ ⎢⎣ ⎡ ⋅ ⋅° Btu or hr ft F W mK . 12 R represents the “short circuit” resistance for heat flow between f 1 T and f 2 T . However, if the fluid temperatures in each leg of the Utube are approximately equal, which occurs at the bottom of the borehole, the resistance 12 R can be neglected in the Δcircuit. The 12 R resistance is often neglected for the entire borehole. This has the effect of decoupling one leg of the Utube from the other, greatly simplifying the system. Figure 14 shows a decoupled borehole system with a circuit diagram defining 1 R and 2 R . GROUT UTUBE GROUND BOREHOLE SYMMETRY LINE RP1 Rf1 Rg1 Rf2 RP2 Rg2 Tf Tb Tb Tf Tf R1 R2 Tb Figure 14 Cross Section of a Borehole with Symmetry Line and the Corresponding Thermal Circuit. In decoupling the borehole, the assumption is made that the grout, pipe, and fluid for each half of the borehole have the same geometry and thermal properties. This 9 assumption means that f 1 f 2 R = R , p1 p2 R = R , and g1 g 2 R = R . Thus 1 R and 2 R , from the circuit in Figure 13, are equal. However, the total resistance of the grout is not typically written in terms of the grout resistance for half of the borehole. The overall grout resistance is instead lumped in one g R term. With these assumptions the borehole resistance circuit shown in Figure 14 can easily be reduced to produce Equation 11. This equation describes the overall borehole resistance. 2 pipe fluid total grout R R R R + = + where, total R = borehole thermal resistance ⎟⎠ ⎞ ⎜⎝ ⎛ W mK or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF grout R = grout thermal resistance ⎟⎠ ⎞ ⎜⎝ ⎛ W mK or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF pipe R = pipe thermal resistance for one tube ⎟⎠ ⎞ ⎜⎝ ⎛ W mK or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF fluid R = fluid thermal resistance for one tube ⎟⎠ ⎞ ⎜⎝ ⎛ W mK or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF (11) The two major contributors to the borehole resistance are the grout and pipe resistance. The fluid resistance contributes typically less than one percent to the overall steady state borehole resistance for turbulent flow. For laminar flow the contribution made by the fluid resistance is much greater and can exceed twenty percent of total R . The pipe resistance can be calculated with Equation 1–2 (Drake and Eckert, 1972). 10 k r r Rpipe 2π ln 1 2 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ = where, pipe R = pipe thermal resistance ⎟⎠ ⎞ ⎜⎝ ⎛ W mK or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF k = Conductivity of the pipe ⎟⎠ ⎞ ⎜⎝ ⎛ mK W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h⋅ ft⋅ F Btu o 2 r = outside diameter (m) or (ft) 1 r = Inside diameter (m) or (ft) (12) The fluid resistance can be calculated using Equation 13 (Drake and Eckert, 1972). r h Rfluid 1 2 1 π = where, fluid R = fluid thermal resistance ⎟⎠ ⎞ ⎜⎝⎛ W mK or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF h = convection coefficient of the fluid ⎟⎠ ⎞ ⎜⎝ ⎛ m K W 2 or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft ⋅ F Btu 2 o 1 r = Utube inside diameter (m) or (ft) (13) The grout resistance can be calculated from the average temperature profile at the borehole wall and the surface of the Utubes with Equation 14 (Hellström, 1991), presuming these temperatures are available. 11 Q R TU tube TBH Wall grout − − − = where, grout R = grout thermal resistance ⎟⎠ ⎞ ⎜⎝ ⎛ W mK or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF Q = Heat flux per unit length of Utube ⎟⎠ ⎞ ⎜⎝ ⎛ m W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft Btu U tube T − = Average temperature at outer surface of Utube (K) or (ºF) BH Wall T − = Average borehole wall temperature (K) or (ºF) (14) The grout resistance is the most complicated component of the borehole resistance. Unlike the pipe and fluid resistance, the apparent grout resistance will change significantly over the first few hours of heat injection or extraction. There are several methods of calculating the grout thermal resistance. The methods that are used to directly calculate the steady state borehole resistance are the Gu and O’Neal (a, 1998) approximate diameter equation, the Paul (1996) method, and the multipole (Bennet and Claesson, 1987) method. Other methods calculate the transient heat transfer between the fluid and surrounding ground, but may be applied to calculate the borehole resistance, include the cylinder source (Ingersoll, 1948) method, and the finite volume method (Yavuzturk, 1999). 1.2.1.1 Line Source Model The line source developed by Kelvin and later solved by Ingersoll and Plass (1948), is the most basic model for calculating heat transfer between a line source and the earth. In this model the borehole geometry is neglected and modeled as a line source or sink of infinite length, surrounded by an infinite homogenous medium. Thus, with 12 respect to modeling a borehole, the line source model neglects the end temperature effects. Equation 15 is the general equation that Ingersoll and Plass (1948) used to model the temperature at any point in an infinite medium from a line source or sink. The medium is assumed to be at a uniform temperature at time zero. ∫ ∞ − Δ = soil x e d k T q β π β β 4 where, t x r soil 4α 2 = (15) (16) ΔT = change in ground temperature at a distance r from the line source (°C) or (°F) q = heat transfer rate per length of line source ⎟⎠ ⎞ ⎜⎝ ⎛ m W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h⋅ ft Btu t = time duration of heat input q (s) r = radius from the line source (m) or (ft) soil α = soil thermal diffusivity ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ s m2 or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ s ft2 soil k = conductivity of the soil ⎟⎠ ⎞ ⎜⎝⎛ mK W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft⋅ F Btu o The integral in Equation 15 can be approximated with Equation 17 for an mth stage of refinement. 13 Σ= ⋅ − = − − − m n n n n n I x x x 1 ! ( ) γ ln( ) ( 1) where, x = Defined by Equation 1–6 γ = 0.5772156649 = Euler’s Constant (17) Equation 17 shows the general form of the line source for the mth stage of refinement. In most references the second stage of refinement is used. This method is only accurate for large times. For a typical borehole this equates to times greater than approximately 10 hours. For small times, less than 10 hours, the GaussLaguerre quadrature approximation, as shown in Equation 18, is given. This approximation uses the fourth order GaussLaguerre quadrature to solve the infinite integral in Equation 15. ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ + + ⋅ + + ⋅ + + ⋅ + = − ⋅ 44 44 43 43 42 42 41 41 ( ) 1 1 1 1 x z w x z w x z w x z I x e x w quad where, 41 w = 0.6031541043 41 z = 0.322547689619 42 w = 0.357418692438 42 z = 1.745761101158 43 w = 0.0388879085150 43 z = 4.536620296921 44 w = 0.00053929470561 44 z = 9.395070912301 (18) The Gauss – Laguerre quadrature approximation shown here should be used for small times, approximately less than 10 hours. Equation 19 is a modification of Equation 15 and shows how to use the line source to model the borehole fluid temperature. Without the borehole resistance ( bh q ⋅ R ) 14 the borehole temperature (T ) would be the temperature at the borehole wall radius and not the fluid temperature. bh ff soil bh soil q R T t I r k T t q + ⋅ + ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ = ⋅ 4π 4α ( ) 2 (19) T = Borehole fluid temperature (°C) or (°F) t = time duration of heat input (s) ff T = far field temperature of the soil (°C) or (°F) q = heat transfer rate per length of line source ⎟⎠ ⎞ ⎜⎝ ⎛ m W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h⋅ ft Btu bh r = radius of the borehole (m) or (ft) bh R = steady state borehole resistance ⎟⎠ ⎞ ⎜⎝ ⎛ W mK or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF soil α = soil thermal diffusivity ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ s m2 or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ s ft2 soil k = conductivity of the soil ⎟⎠ ⎞ ⎜⎝ ⎛ mK W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft⋅ F Btu o Equation 19 differs from Equation 1.5 in that it uses the steady state borehole resistance to model the heat transfer from the borehole wall to the fluid; the line source model is used to model the heat transfer between the borehole wall and the far field. This usage of the line source requires that the steady state borehole resistance is known. Since the steady state resistance is used the line source will have error for short times before the borehole reaches steady state resistance. For most boreholes the error in the steady state borehole resistance is negligible at 2 hours. The line source model is very easy to use and requires relatively few calculations compared to other methods. However, the drawback to this model is that the borehole 15 internal geometry, thermal properties, and the mass of the fluid are not modeled. The resulting inaccuracies will be examined in Chapter 4. 1.2.1.2 GuO’Neal Equivalent Diameter Model The Gu and O’Neal (1998 a) equivalent diameter method is a very simple method of calculating the steady state borehole thermal resistance. It yields a steady state borehole resistance value that is adequate for most simple calculations. This method is represented by an algebraic equation for combining the Utube fluid into one circular region inside the center of the borehole such that the resistance between the equivalent diameter and borehole wall is equal to the steady state borehole resistance of the grout. Equation 1–10 is used to calculate the equivalent diameter. As can be seen the equivalent diameter is based solely on the diameter of the Utube and the center to center distance between the two legs. eq s D = 2D⋅ L s BH D ≤ L ≤ r where, eq D = Equivalent diameter (m) or (ft) BH r = radius of the borehole (m) or (ft) D = diameter of the Utube (m) or (ft) s L = center to center distance between the two legs (m) or (ft) (110) Figure 15 shows three actual configurations and their equivalent diameters. “d” shows the equivalent diameter for configuration “a”; “e” for “b”; and “f” for “c”. 16 (a) (b) (c) (d) (e) (f) Figure 15 Actual Geometry vs Equivalent Diameter Approximation To calculate the grout resistance, Equation 111, which is the general equation for radial heat conduction through a cylinder, should be employed. grout eq bh grout k D D R 2π ln ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = where, grout R = grout thermal resistance ⎟⎠ ⎞ ⎜⎝ ⎛ W mK or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF grout k = conductivity of the grout ⎟⎠ ⎞ ⎜⎝ ⎛ mK W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h⋅ ft⋅ F Btu o bh D = diameter of the borehole (m) or (ft) eq D = equivalent diameter using GuO’Neal’s method (m) or (ft) (111) The steady state resistance of the grout can be used in Equation 1.1 to calculate the overall borehole resistance. 17 1.2.1.3 Paul Model An experimentally and analytically based method for calculating steady state borehole resistance was developed at South Dakota State University by Paul (1996). The Paul method for calculating the steady state borehole resistance was created using both experimental data and a two dimensional finite element program for modeling a borehole cross section. Several different borehole parameters were modeled such as shank spacing, borehole diameter, Utube diameter, grout conductivity, and soil conductivity. The test apparatus used a single layer thick coil of wire wrapped around each side of the Utube to form an electrical resistance heater. This provided a uniform, constant flux heat input for the system. A real borehole will not have uniform flux at the pipe wall. Heat was input until steady state temperature conditions at the borehole wall radius and along the circumference of the Utube were reached. The borehole resistance was then calculated from the temperatures and the flux. A two dimensional finite element model was created using ANSYS, a UNIX based software package, for the purpose of extending the range of borehole diameters and pipe sizes that the steady state borehole resistance could be solved for. The ANSYS cases could be run much faster than the experimental apparatus; this allowed for more cases to be run. Experimental results from the test apparatus and the ANSYS model were compared for validation purposes. From the results, shape factor correlations were created to model the complex geometry of the borehole. Equation 112 is the resulting shape factor equation for calculating the steady state grout resistance. 18 K S R grout grout ⋅ = 1 , 1 0 β β ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ = ⋅ U−tube b d S d where, grout R = Equivalent diameter (m) or (ft) grout K = Conductivity of the grout ⎟⎠ ⎞ ⎜⎝ ⎛ mk W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft⋅ F Btu o S = Shape factor (dimensionless) 0 β and 1 β = Curve fit coefficients (dimensionless) b d = Diameter of the borehole (m) or (ft) U tube d − = Outside diameter of the Utube (m) or (ft) (112) Equation fit coefficients are given by Paul (1996) for four different shank spacings; A0, A1, B and C. The shank spacings are described in Figure 16. Index Spacing Condition A0 S1 = 0 A1 S1=.123 in (.3124 cm) B S1 = S2 C S2 = 0.118 in (.300 cm) S1 S2 Figure 16 Types of shank spacing used in the Paul borehole resistance approximation. The 0 β and 1 β values for the shank spacing described in Figure 16 are given in Table 1 6. Table 11 Paul Curve Fit Parameters used to Calculate the Steady State Grout Resistance A0 A1 B C 0 β 14.450872 20.100377 17.44268 21.90587 1 β 0.8176 0.94467 0.605154 0.3796 R 0.997096 0.992558 0.999673 0.9698754 19 The R value indicates the accuracy of the curve with respect to the experimental or ANSYS model. An R value of 1 indicates a perfect fit. The grout resistance found using this method should be applied within Equation 11 to determine the overall borehole resistance. 1.2.1.4 Cylinder Source Model The cylinder source model, created by Ingersoll and Plass (1948), uses an infinitely long cylinder inside an infinite medium with constant properties and solves the analytical solution of the 2D heat conduction equation. The cylinder source solution for the gfunction and temperature change at the borehole wall can be calculated with Equations 113, 14, and 15. ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ Δ = ⋅ o o soil r G F r k T q , , r 2 t F soil o ⋅ = α ( ) ( ) ( ) ( ( ) ( )) β β β β β β β β π β d J Y r Y J Y r r J r e r G F r F o o o o o ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ + ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ − ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ − = ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ∫ ∞ − ⋅ 2 2 1 2 1 0 1 1 0 0 2 , 1 1 2 where, (113) (114) (115) ΔT = temperature difference between the steady state temperature of the ground and the temperature at the borehole wall (ºC) or (ºF) q = heat flux per unit length of the borehole ⎟⎠ ⎞ ⎜⎝ ⎛ m W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft Btu o F = Fourier number (dimensionless) r = inner cylinder radius (equivalent Utube radius) (m) or (ft) o r = outer cylinder radius (borehole radius) (m) or (ft) 0 J , 1 J , 0 Y , 1 Y = Bessel functions of the zero and first orders t = time (s) 20 The variable “r” is the location at which a temperature is desired from the cylinder source located at o r . G( o F , r / o r ) is a function of time and distance only. To apply the cylinder source equation for modeling the fluid temperature within a borehole, Equation 116 can be used, setting r equal to the equivalent Utube radius and o r equal to the borehole radius. bh ff o o soil q R T r G F r k T t q + ⋅ + ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ( ) = ⋅ , where, (116) T(t) = borehole fluid temperature (ºC) or (ºF) q = heat flux per unit length of the borehole ⎟⎠ ⎞ ⎜⎝ ⎛ m W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft Btu o F = Fourier number (dimensionless) r = inner radius (equivalent Utube radius) (m) or (ft) o r = outer cylinder radius (borehole radius) (m) or (ft) ff T = far field temperature of the soil (ºC) or (ºF) t = time (s) bh R = steady state borehole resistance ⎟⎠ ⎞ ⎜⎝ ⎛ W mK or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF soil k = soil conductivity ⎟⎠ ⎞ ⎜⎝ ⎛ mK W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF The cylinder source can be used to model the steady state borehole resistance using Equation 11. The Utube and fluid resistances can be calculated as shown in section 1.2.1 in Equations 12 and 13. 21 A comparison of borehole resistance calculation methods, including the cylinder source method, is shown in chapter 2. Also a correction factor was created that increases the accuracy of the cylinder source greatly for boreholes with a large shank spacing. As will be shown in chapter 2, even without the correction factor, the cylinder source with superposition technique is a reasonably accurate method for calculating the grout resistance. 1.2.1.5 Multipole The multipole (Bennet, et al. 1987) method is used to model conductive heat flow in and between pipes of differing radius. In the multipole model, the tubes are located inside a homogenous circular region that is inside another homogenous circular region. The multipole method is not constrained to calculating the steady state borehole resistance for a borehole with only one Utube. Furthermore, the tubes do not need to be symmetrical about any axis. This is advantageous since some boreholes have two Utubes. The model is also able to calculate borehole resistance for Utubes that are not equidistant from the center of the borehole. To show the capabilities of the model Figure 17 has been created showing an asymmetric borehole with three pipes. The pipes have temperatures, f 1 T , f 2 T , and f 3 T . 22 Soil Grout Pipe Fluid Tf1 Tf2 Tf3 Ks 1 Kg b 2 b b3 x rb rs 1 r 2 r r 3 bc Tc Figure 17 Example of a 2D System for the Multipole Method. The inner circular region represents the grout and the outer region represents the soil for the borehole system. For calculating borehole resistance, b r can be set to 100 m (328 ft). The inputs to the multipole method are shown in Table 12. 23 Table 12 Variable Input List for the Multipole Method. g K Thermal conductivity in the inner region ⎟⎠ ⎞ ⎜⎝ ⎛ mk W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft⋅ F Btu o s K Thermal conductivity in the outer region ⎟⎠ ⎞ ⎜⎝ ⎛ mk W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft⋅ F Btu o N Number of pipes J Order of multipole b r Radius of the outer region (m) or (ft) s r Radius of the inner region (m) or (ft) c β Thermal resistance coefficient at the outer circle (nondimensional) c T Temperature of the outer region (K) or (ºF) The following are input for each pipe indexed by i i i x , y Location of each innermost pipe (m) or (ft) i r Radius of each pipe (m) or (ft) i β Thermal resistance coefficient for each pipe (nondimensional) fi T Fluid temperature (K) or (ºF) In Table 12 the nondimensional variable i β is used to input the pipe thermal resistances. This is shown in Equation 117. 24 β = 2πRk where, R = Thermal resistance of the pipe ⎟⎠ ⎞ ⎜⎝ ⎛ W Km or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF k = Grout conductivity ⎟⎠ ⎞ ⎜⎝ ⎛ mk W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft⋅ F Btu o β = resistance coefficient (Nondimensional) (117) The general equation that the multipole method solves is the steady state twodimensional heat conduction equation, Equation 118: 0 2 2 2 2 = ∂ ∂ + ∂ ∂ y T x T where, T = Temperature (°C) or (ºF) x = Distance in the x direction (m) or (ft) y = Distance in the y direction (m) or (ft) (118) When solving the differential equation several assumptions are made. The temperature is constant inside of the pipe walls and the fluid convective resistance. The temperature around the outer region is constant. The system is at steady state. Multipoles were created using the line source model. They are called multipoles because for each line source there is a line sink at a mirror point. This can be seen in Equation 119, the temperature equation for a zeroth order multipole, where the line sink n q is at ( , ) n n x y in the first term and the mirror sink of strength n σ ⋅ q is located at the mirror point ( 2 / 2 , 2 / 2 ) n b n n b n x r r y r r in the second term. A zeroth order multipole has one source and one sink. 25 0 ≤ r ≤ rb : (119) ⎪ ⎪ ⎭ ⎪ ⎪ ⎬ ⎫ ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ − + − + ⋅ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − + − = 2 2 2 2 2 2 2 2 ( ) ( ) ln / ( ) ( ) ln 2 ( , ) n b n n b n b n n n b b n r y y r r x x r r r x x y y r k T x y q σ π . k k k k b b + − σ = , −1≤σ ≤1 where, T(x, y) = Temperature (K) or (ºF) n q = Heat flux per unit length ⎟⎠ ⎞ ⎜⎝ ⎛ m W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft Btu b k = Conductivity of the inner region ⎟⎠ ⎞ ⎜⎝ ⎛ mk W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft⋅ F Btu o k = Conductivity of the outer region ⎟⎠⎞ ⎜⎝ ⎛ mk W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft⋅ F Btu o n x = Location of the line source in the x direction (m) or (ft) n y = Location of the line source in the y direction (m) or (ft) b r = Radius of the inner region (m) or (ft) n r = Radius of the line source (m) or (ft) To give a graphical perspective on the location of the source and sink Figure 1.8 is given. As can be seen in Figure 1.8 the sink lies on the same radius line as the source. 26 Sink Source y x Grout Soil Radial line from borehole center Borehole Pipe Figure 18 Source and Sink Location for a Single Pipe Equation 119 is the zeroth order multipole equation; to produce a more accurate solution more sources and sinks can be utilized for each pipe. To do this requires a simplification of Equation 119 by use of polar notation. Writing Equation 119 in polar notation yields Equation 120 where ( ) 0 Re n W is the real component of the zero order multipole. 27 ( ) 0 Re 2 ( , ) n b n W k T x y = q ⋅ π where, T(x, y) = Temperature (K) or (ºF) n q = Heat flux per unit length ⎟⎠ ⎞ ⎜⎝ ⎛ m W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft Btu b k = Conductivity of the inner region ⎟⎠ ⎞ ⎜⎝ ⎛ mk W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft⋅ F Btu o ( ) 0 Re n W = Real component of the zero order multipole (120) For higher order multipoles, derivatives are taken of the n0 W as shown in Equation 121. ( ) ( 1)! 1 j no n j nj W j z W ∂ ∂ ⋅ − = where, nj W = jth order multipole of nth line source j = Order of multipole n z = Location of pipe n in polar coordinates no W = zero order multipole (121) Both the real and imaginary components of nj W satisfy the continuous radial flux boundary condition at b r = r . The constant temperature condition of each of the pipes and the outer radius s r is satisfied using a Fourier series expansion. Using this method the temperature for any point within s r can be found. Equation 122 and 23 show the 28 general solution for the temperature inside the outer cylinder radius for orders of multipoles greater than zero. For a borehole this becomes the borehole wall temperature. ⎥⎦ ⎤ ⎢⎣ ⎡ = + Σ ⋅ +ΣΣ ⋅ ⋅ +Σ ⋅ − ⋅ = = j cj j nj cj c N n j j nj pn N n o n no T T P W P r W P r W 1 1 Re b n n k P q ⋅ = 2π where, (122) (123) T = Temperature at the borehole wall (K) or (ºF) o T = Far field temperature (K) or (ºF) n = counter variable j = Order of multipole N = number of pipe no W = zero order multipole cj W = Multipole of the outer cylinder nj W = jth order Multipole of nth line source pn r = Radius of the innermost pipes (m) or (ft) c r = Radius of the cylinder encircling the innermost pipes (m) or (ft) The final equations, shown in Chapter 8 (Bennet, et al. 1987), are an elaboration of Equation 122. In the paper three equations are presented which must be solved iteratively. In Chapter 11 (Bennet, et al. 1987) Fortran 77 code is conveniently given which solves the equations. 29 The multipole method can produce highly accurate steady state temperature profiles of a borehole and soil. Figure 19 shows the two dimensional steady state temperature inside and around a typical borehole. For this figure the borehole fluid temperature was set to 10 °C (18 °F) above a zero far field temperature and a tenth order multipole was used. This method requires a constant temperature boundary condition inside the fluid convective resistance inside the Utubes. This is a reasonable assumption if the fluid in the Utubes is in the turbulent flow regime. Figure 19 Steady State Temperature Field for a Borehole Heat Exchanger 30 A tenth order multipole produces four or five digits of accuracy. Since the multipole method is very fast on computers above 400 megahertz, a tenth order multipole should be used for most problems. As can be seen in Figure 19, the difference in grout and soil conductivity creates a slope discontinuity at the borehole radius, where the inner circular region representing the grout meets the outer region representing the soil. A typical steady state temperature profile at the borehole wall is shown in Figure 110. Steady State Borehole Wall Temperature 0 1 2 3 4 5 6 7 8 0 50 100 150 200 250 300 350 Distance Along the Circumference (Degrees) Temperature (c) Figure 110 Temperature Change Around the Borehole Circumference To calculate the steady state borehole resistance a temperature is specified for each leg of the Utube. The multipole method is then used to calculate a heat flux out of each Utube and the temperature around the circumference of the borehole radius should be averaged. Since the multipole program solves for temperatures very quickly, using 360 points at the borehole radius is feasible and will produce a high degree of accuracy. 31 Once an average temperature at the borehole radius is found and it can then be used with the flux in Equation 124 to find the steady state borehole resistance. U tube1 U tube2 fluid bh BH q q T T R − − + − = where, BH R = borehole thermal resistance ⎟⎠ ⎞ ⎜⎝ ⎛ W mK or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF fluid T = average temperature at the Utube pipe wall (K) or (ºF) bh T = average temperature at the borehole radius (K) or (ºF) U tube1 q − = heat flux from Utube leg one ⎟⎠ ⎞ ⎜⎝ ⎛ m W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft Btu U tube2 q − = heat flux from Utube leg two ⎟⎠ ⎞ ⎜⎝ ⎛ m W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft Btu (124) 1.2.2 Transient Modeling of Borehole Heat Exchangers Transient heat transfer in boreholes occurs when the heat flux entering the borehole through the fluid does not equal the heat flux leaving the borehole via the borehole wall. Borehole transients have a significant effect on the borehole fluid temperature response after any change in heat extraction or rejection rate. For a step change in the heat flux, the time for which transient effects are significant is determined primarily by the grout thermal conductivity and the borehole geometry such as the shank spacing and radius of the borehole. In general, a small grout thermal conductivity or a large borehole radius will lengthen the transient region. In most actual systems, the heat flux applied to a borehole through the circulating fluid changes continuously throughout operation. Thus borehole transients need to be modeled not only at the beginning of a simulation but throughout the simulation. 32 Several analytical two dimensional models exist and have been used to model boreholes such as the line source (Ingersoll and Plass, 1948) and cylinder source (Ingersoll, 1948). These methods have very limited capability for modeling the internal borehole transients especially for transient heat pulse changes in the first 10 hours. This section, describes the buried electrical cable (Carslaw and Jaeger, 1947) analytical model, and presents the General Elliptical Multiblock Solver (GEMS2D) and Yavuzturk’s (1999) pie sector finite volume method programs. The application of the buried electrical cable model to borehole heat exchangers is covered in Chapter 3. This section also covers Eskilson’s (1987) three dimensional model of boreholes and its coupling with the two dimensional analytical or finite volume methods via borehole resistance. 1.2.2.1. Buried Electrical Cable Model The buried electrical cable (BEC) model (Carslaw and Jaeger 1947) is an analytical model used to describe the heat flow out of a cable buried in the ground. An electrical cable consists of three main parts, a metal core surrounded by insulation and then an outer protective sheath. A diagram of a buried electrical cable is shown in Figure 111 along with a circuit diagram of the system. Implicit in this method are the assumptions that the core and sheath thermal capacities have finite thermal capacities but are perfect conductors and that the insulation has negligible heat capacity, but a fixed thermal resistance. The most significant difference between the buried electrical cable model and other analytical models such as the line source and cylinder source is that this model incorporates the thermal capacity of the sheath and core in calculating the temperature profile of the core. As seen in the circuit diagram in Figure 111 the core and sheath 33 thermal capacities are represented as 1 S and 2 S . The insulation resistance is represented as s R . In the circuit a heat flux can be applied creating a temperature differential between the core and soil. The heat flux is absorbed by the capacitances 1 S and 2 S . INSULATION SHEATH GROUND CORE Rs Tsoil Tcore Tsoil, ff Tsoil, ff S1 S2 Figure 111 Diagram of a Buried Electrical Cable and Circuit The analytical equations, given as Equations 125 and 126, for this system are more complicated and require more computational time than the cylinder source and the line source equations. 34 G(t) k T q soil Δ = (125) du u G t a a e a t u ∫∞ ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ − ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ⋅ Δ − = 3 0 3 2 2 2 1 2 2 ( ) 2 1 α π (126) 1 2 1 2 S r C a b p π ρ = , 2 2 2 2 S r C a b p π ρ = , s soil h = 2π ⋅ R k [[ ( ) ( ) ] [ ( ) ( ) ]2 ] 1 2 0 2 1 2 1 2 2 1 2 0 2 1 2 1 2 Δ = u a + a − hu J (u) − a a − hu J (u) + u a + a − hu Y (u) − a a − hu Y (u) where, t = time (s) b r = outer radius of the sheath (m) or (ft) soil k = Conduc1tivity of the soil ⎟⎠ ⎞ ⎜⎝ ⎛ mK W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft⋅ F Btu o ρ = density of the soil ⎟⎠ ⎞ ⎜⎝ ⎛ m3 kg or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ft3 lbm p C = specific heat of the soil ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ kgK J or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ kg⋅ F Btu o s R = insulation thermal resistance ⎟⎠ ⎞ ⎜⎝ ⎛ W Km or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF 1 S = core thermal capacity ⎟⎠ ⎞ ⎜⎝ ⎛ mK J or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ft⋅ F Btu o 2 S = sheath thermal capacity ⎟⎠ ⎞ ⎜⎝ ⎛ mK J or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ft⋅ F Btu o soil α = soil thermal diffusivity ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ s m2 or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ s ft2 0 J , 1 J , 0 Y , 1 Y = Y and J type Bessel functions of zero and first orders 35 Equation 126 will produce a buried electrical cable gfunction for a particular time. It should be noted that this is the only analytical model presented here which takes into account the thermal mass of the heat generation medium, which in the case of a buried electrical cable is the core. However, there is potential for this model to be modified to model a borehole and account for the fluid mass inside a borehole. The application of the BEC analytical equation in modeling a borehole system is discussed in detail in Chapter 3. 1.2.2.2. Gfunction Model: Long Time Step The long time step (LTS) gfunction (Eskilson 1987) represents the nondimensionalized borehole response for times when threedimensional effects such as borehole to borehole interaction, surface and bottom end effects influence the borehole fluid temperature response. Gfunctions are plotted against the natural log of scaled time where the scaling factor is dependent on the depth of the borehole and the soil thermal diffusivity. As developed by Eskilson (1987), the borehole transient resistance or gfunction can be nondimensionalized with respect to the soil and scaled with respect to the steady state borehole resistance to form a gfunction. Both long and short time step gfunctions can be produced using Equation 127 (Eskilson 1987). In Equation 127 the borehole T term represents the timevarying average temperature at the borehole wall and must be calculated with a numerical or analytical procedure. ground T is the far field temperature and usually remains constant. This gfunction represents the nondimensionalized resistance between the ground and borehole wall. Equation 128 includes the borehole resistance term. 36 ( borehole ground ) b soil s T T Q k H r t g t = − ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ 2π , ( ) borehole ground soil BH b soil s T T k R Q k H r t g t π π 2 2 , + − = ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ where, g = gfunction value (dimensionless) Q = flux per unit length ⎟⎠ ⎞ ⎜⎝ ⎛ m W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft Btu soil k = thermal conductivity of the soil ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ft⋅ F Btu o borehole T = average temperature at the borehole wall (°C) or (°F) ground T = far field temperature of the ground (°C) or (°F) (127) (128) As can be seen from Equation 128, the gfunction has two major parameters s t / t and r H b / . For a specific borehole configuration, the first parameter is the major contributor to the gfunction and the second is a factor that corrects the gfunction according to the borehole radius ( b r ) to depth (H). The r H b / correction factor is relatively minor since it changes the gvalues on the order of one percent or less. The main parameter, in Equation 127, that requires significant calculation time is the average temperature at the borehole wall radius ( borehole T ). Gfunctions are plotted against the natural log of nondimensionalized time. The term s t is called the time scale factor, and can be calculated using Equation 129 (Eskilson 1987). 37 soil s t H 9α 2 = where, s t = time scale factor (s) H = depth of the borehole (m) or (ft) soil α = soil thermal diffusivity ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ s m2 or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ s ft2 (129) The gfunction defined by Equation 127 represents the ground thermal resistance and not the borehole resistance. This is beneficial because it allows a single long time step gfunction to be useful for any borehole geometry and soil conductivity as discussed by Eskilson (1987). The gfunction in Equation 128 is only valid for the specific borehole for which the borehole resistance was calculated. An example of a combined long and short time step gfunction for a single borehole system is shown in Figure 112 and 113. In these figures, the gfunction is plotted against log scale time. Figure 112 was created with Equation 127 and Figure 1 13 with Equation 128. Short and Long Time Step Gfunction for a Single Borehole System Without Borehole Resistance 3 2 1 0 1 2 3 4 5 6 7 16 14 12 10 8 6 4 2 0 2 4 ln(t/ts) gfunction Figure 112 Short and long time step gfunction without borehole resistance Short and Long Time Step Gfunction for a Single Borehole System with Borehole Resistance 0 1 2 3 4 5 6 7 8 9 10 16 14 12 10 8 6 4 2 0 2 4 ln(t/ts) gfunction Figure 113 Short and long time step gfunctions with borehole resistance 38 As can be seen in Figure 113 the gfunction approaches zero at small times. This indicates that as time approaches zero the resistance asymptotically approaches zero due to the steady state temperature profile of the soil. Looking at Figure 112 might give the impression, however, that for small times the resistance is negative but this is an illusion created by subtracting the borehole resistance. When the borehole resistance is added back in, as shown in Figure 113, the resistance approaches zero at short times. At large times, Gfunctions will plateau. This occurs because of borehole end effects. Using Equation 130 an average fluid temperature can be calculated if the gfunction is known. ground b soil s borehole T H r t g t k T Q + ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ = , 2π where, borehole T = average temperature at the borehole wall radius (ºC) or (ºF) ground T = far field temperature of the ground (ºC) or (ºF) g = gfunction value (dimensionless) Q = flux per unit length ⎟⎠ ⎞ ⎜⎝ ⎛ m W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ hr ⋅ ft Btu soil k = thermal conductivity of the soil ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ft⋅ F Btu o (130) There are no published analytical solutions that approximate the gfunctions for multiple borehole systems. This is due to multiple borehole systems dependence on not only the depth of the borehole, but also on the distance between each borehole. The interaction between boreholes is difficult or impossible to analytically model. With the 39 boreholes dissipating different unknown amounts of heat, it is difficult to analytically model the average borehole wall temperature of the system. However, it can be resolved using superposition. In order to create long time step gfunctions for multiple borehole systems Eskilson (1987) created a Fortran 77 program which uses a variable mesh finite difference method with cylindrically symmetric coordinates. The program is described in detail for a single borehole in Eskilson (1987). The program has the ability to input a constant heat flux per unit length of borehole and calculates the resulting temperature at the borehole radius ( borehole T ) for various times. The borehole wall temperature ( borehole T ) is then used in Equation 127 to calculate the LTS gfunction. An example of the type of variable mesh grid that Eskilson used is shown in Figure 114 Borehole (i1, j) (i, j) (i, j1) (i+1, j) (i, j+1) Figure 114 Twodimensional radialaxial mesh for a heat extraction borehole in the ground (Eskilson, 1987) In Figure 114 each cell is a rectangular cross section of a ring. Temperatures are calculated at the center of each cell; however, logarithmic interpolation can be used to 40 find the temperature at other points in the soil. Eskilson gives a detailed analysis of appropriate mesh sizes in the radial and axial direction since the mesh determines the accuracy of the solution. In general, small cells provide good numerical accuracy and the temperatures are valid for smaller times. However, small cells create longer computational times for a computer. Eskilson suggests using no smaller cells than necessary for a particular problem. Although computer technology has improved since 1987, mesh size is still important for large simulations. In the examples Eskilson used, the upper part of the borehole is thermally insulated to a depth of 5 meters (16.4 ft) and the overall borehole depth is 115 m (377 ft). Mesh comparisons for short times of 25 years and long times of 237 and 947 years were conducted. In the end, heuristics were created for determining the appropriate mesh size in the radial and axial direction. To solve the thermal performance for a system with multiple boreholes, Eskilson (1987) used the superposition technique. Two different examples are given (Eskilson, 1987): a 4x4 borehole configuration, with 10 m (33 ft) spacing, and a 12x10 borehole configuration, with 4 m (13 ft) spacing, with the simplest type of loading condition where the heat flux at the borehole wall is constant per unit length of the borehole. To validate the accuracy of the program the line source was used in conjunction with superposition. Eskilson’s program was used to produce the LTS gfunction curves shown in Figure 115. It can intuitively be determined that a tighter borehole field will produce more overall resistance and as the boreholes are spaced farther apart, all multiple borehole systems will approach the single borehole case. This can be seen in Figure 1 15. 41 Long Time Step Gfunction for Different Borehole Spacing 0 10 20 30 40 50 60 70 80 90 100 9 8 7 6 5 4 3 2 1 0 1 2 3 ln(t/ts) gfunction 1.5 m (5ft) 3.0 m (10 ft) 4.6 m (15 ft) 6.1 m (20 ft) 9.1 m (30 ft) 24 m (80 ft) Figure 115 Long Time Step gfunction for a 64 Borehole System in an 8x8 Configuration with Varying Borehole Spacing Since none of the internal properties of the borehole are significant, the long time step gfunction might initially seem simple to solve. Because of three dimensional effects, the gfunction for multiple borehole systems is deceptively complicated to solve. 1.2.2.3. Gfunction Model: Short Time Step The short time step (STS) gfunction describes the transients that occur within the borehole before the borehole reaches steady state conditions. For this transient region, the borehole is modeled as having infinite length since surface and bottom end effects can be neglected. The STS gfunction can be approximated using the line source or the cylinder source as described in section 1.2.1.1 and 1.2.1.4, respectively, to calculate a fluid temperature profile versus time which could be input into Equation 131 to yield a 42 gfunction. Equation 131 calculates a gfunction in terms of the fluid temperature and the steady state borehole resistance. Unlike the long time step gfunction, this gfunction is only valid for the specific borehole internal geometry (shank spacing, borehole radius… etc.) and conductivities that the fluid T and BH R was generated with. ( fluid BH ground ) b soil s T R Q T Q k H r t g t = − − ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ( ) 2 , π where, g = gfunction value (dimensionless) Q = flux per unit length ⎟⎠ ⎞ ⎜⎝ ⎛ m W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft Btu soil k = thermal conductivity of the soil ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ft⋅ F Btu o ground T = far field temperature of the ground (°C) or (°F) fluid T = average temperature of the circulating fluid (°C) or (°F) BH R = borehole resistance ⎟⎠ ⎞ ⎜⎝ ⎛ W mK or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF (131) Another method which is capable of generating greater accuracy than analytical methods is the finite volume model (Patankar, 1980). Two programs that have implemented this model will be discussed. The first is called the General Elliptical Multiblock Solver (GEMS2D) and was developed by Rees (2001). Applications of GEMS2D are reported by Spitler, et al. (1999) and Rees, et al. (2002). This program solves the general convection diffusion equation using a boundary fitted grid. GEMS2D is capable of solving both steady state and transient problems. Boundary fitted grids enable GEMS2D to be applied in solving heat transfer problems with complex geometries such as Utubes within a borehole. Figure 116 shows a GEMS2D boundary fitted grid for half of a borehole, since the geometry is symmetrical. 43 Figure 116 Grid for a cross section of a borehole A complicated grid such as that shown in Figure 116 will require several different blocks to be created and then connected together. Each block is composed of many cells. The cells in each block can then be assigned properties such as conductivity and heat capacity. Different cells within a single block can be assigned different properties. A detailed description of how the GEMS2D program was applied to borehole heat conduction with fluid mass is given in Section 4.1. GEMS2D is capable of calculating the steady state borehole resistance and the transient temperature profile at the borehole wall. The gfunction can be calculated from the average borehole wall temperature ( borehole T ) using Equation 131. GEMS2D is written in the Fortran 90/95 language. A grid generation tool was also written to automate the creation of grids for the GEMS2D simulator. Using a text input file, the grid for convectiondiffusion heat transfer problems can be created with the grid generation tool. In the text file, blocks and boundaries are created and thermal properties of each block are specified. After the grid is created, GEMS2D can then be 44 used to simulate the system. The GEMS2D outputs are given in an output text file. Temperatures at each node at each time increment can be given for transient simulations. The second program was developed by Yavuzturk, et al. (1999) and also uses the two dimensional finite volume model, but with a polar grid. It is specifically developed for modeling the heat flow out of a Utube borehole heat exchanger. Like GEMS2D, Yavuzturk’s program is able to model both the transient and steady state solutions for the temperature field within and around a borehole. To model the geometry of a borehole heat exchanger Yavuzturk uses an approximation for the borehole Utube geometry called the pie sector approximation (Yavuzturk and Spitler, 1999). The piesector approximates the cross section of the Utubes via two “pieshaped” wedges. Figure 117 shows the grid that Yavuzturk’s program creates. Only half the borehole is shown since the system is symmetrical. The pie shaped wedge shown in this figure is representative of one leg of the Utube. It is shown bolded, in the figure, while the actual circular Utube geometry is shown for both legs without bolding. 45 Figure 117 Grid for PieSector Approximation (Yavuzturk, 1999) The grid resolution and pie sector approximation for the Utube geometry is determined by an automated parametric grid generation algorithm and is a function of the borehole and Utube pipe geometry. The algorithm matches the inside perimeter of the circular pipe to the inside perimeter of the pie sector and also creates identical heat flux and resistance conditions near the pipe wall between the circular pipe and the pie sector approximation. The fluid resistance is approximated by adjusting the thermal conductivity of the Utube pipe wall. The total radius of the grid is 3.6 m or (12 ft) so that longer simulation times can be conducted. At this radius, the boundary condition is set to a constant far field temperature. The model was primarily written in the Fortran 77 programming language. Inputs to the model are simple since grid generation is automated. The inputs include shank spacing, Utube diameter, borehole diameter, convection coefficient, the volumetric heat 46 capacity of the soil, grout and Utubes, and the conductivity of the soil, grout and Utubes. These inputs are provided in a text file. The outputs are given in an output text file. Since the Utube geometry is not modeled as accurately as in GEMS2D, the piesector approximation will generally be less accurate than GEMS2D. The benefit of Yavuzturk’s program is that it requires approximately half the simulation time as GEMS2D. The resulting simulation model, discussed in detail in Yavuzturk and Spitler (1999), uses Eskilson’s LTS gfunctions simulation methodology with Yavuzturk’s STS gfunctions simulation methodology. The model has been incorporated into a commercially available GLHE design tool called GLHEPRO (Spitler, 2000). GLHEPRO Version 3 is discussed in detail in section 1.2.3. The simulation model has also been implemented and proved useful in several studies (Yavuzturk and Spitler, 2000, Ramamoorthy, et al. 2001, and Chiasson, 1999). Hybrid GSHP systems use other heat rejection equipment, such as cooling towers, fluid coolers, shallow ponds (Chiasson, et al. 2000; Ramamoorthy, et al. 2001) or pavement heating systems (Chiasson, et al. 2000). For example, several operating and control strategies of a cooling tower hybrid GSHP system are discussed in Yavuzturk and Spitler (2000), and are compared with hourly simulations performed in TRNSYS (SEL 1997). A goal of these studies is to show that hybrid GSHP systems can reduce the size of the GLHE system which in turn can reduce the first cost of the system and the necessary land area. 47 Spitler, et al. (2000) which gives a summary of research and developments in grounds source heat pump systems, design, modeling, and applications for commercial and institutional buildings. Presented in this paper are design methodologies for determining hourly and minutely responses for GLHE designs. 1.2.3 GLHEPRO Version 3 Design Tool GLHEPRO Version 3 (Spitler, 2000) combines a Microsoft windows graphical user interface and a ground loop heat exchanger simulation. The software package developed by Marshall and Spitler is based on the methods developed by Eskilson (1987) at the University of Lund, Sweden. GLHEPRO Version 3 has a library of heat pump performance curves and has the capability of adding user defined heat pump performance curves. GLHEPRO Version 3 also has the flexibility of using SI or English units. GLHEPRO Version 3 uses the long time step gfunctions developed by Eskilson (1987). As discussed earlier, Eskilson (1987) created a 2 dimensional finite difference program that omits the internal borehole properties such as shank spacing, grout and Utubes properties. This was done by assuming a steadystate heat transfer process inside the borehole and modeling the transient process outside the borehole using a finite difference technique, so that the temperature at the borehole wall ( borehole T ) was found. Equation 131 can then be used to create the gfunction. GLHEPRO Version 3 uses data from Eskilson’s finite difference program to model the long time step gfunction for over 250 different borehole configurations. The fluid temperature is found by using the temperature at the borehole wall in conjunction with the borehole steady state resistance. The method used in GLHEPRO Version 3 to calculate the short time step response is the 48 line source (Ingersoll, 1948). The Paul (1996) method was used to calculate the borehole resistance for a specific borehole geometry. GLHEPRO Version 3 received inputs for peak and monthly, heating and cooling loads along with borehole internal geometric and thermal properties and borehole configuration. GLHEPRO Version 3 has the capability of outputting the maximum and minimum monthly fluid temperature entering the heat pump and the energy consumption of the system. GLHEPRO Version 3 also has a sizing mode which requires maximum and minimum limits for the entering fluid temperature to the heat pump. The sizing program will find the minimum depth required for a specific borehole configuration to be within the maximum and minimum user defined fluid temperature limits. 1.3 OBJECTIVES The primary objective is to develop and implement a method whereby engineers can more accurately model peakloaddominant systems without time consuming numerical modeling. Implicit within this main objective are the following specific objectives: 1. Determine an appropriate method for calculating the steady state borehole resistance and implement it in GLHEPRO. 2. Enhance shorttimestep (STS) GLHE simulation methodology to account for thermal mass of the fluid to yield more accurate designs via simulations. 3. Develop an automated method for producing the combined short and long time step gfunction. 4. Evaluate the impact of the more accurate gfunction calculation methodologies on the simulation of GLHE systems. 49 5. Evaluate the impact of the more accurate gfunction methodologies on the simulation of GSHP. 50 2 COMPARISON OF BOREHOLE RESISTANCE CALCULATION METHODS A small change in the steady state borehole resistance has a significant impact on the borehole fluid temperature profile. Since the short time step (STS) gfunction is derived directly from the borehole fluid temperature it includes the borehole resistance. In order to be consistent with the long time step (LTS) gfunction, it is necessary to adjust the STS gfunction to subtract the nondimensional temperature rise due to the borehole resistance. This, in turn, requires accurate knowledge of the borehole resistance. Furthermore, as discussed in Chapter 3, the borehole fluid thermal mass (BFTM) model, which is developed in this thesis for calculating the STS gfunction, requires the steady state borehole resistance to be known. Thus, there was a need for comparing different borehole resistance calculation methods for the purpose of choosing one for the BFTM model. This chapter provides a comparison between different methods for calculating the steady state borehole thermal resistance. Both numerical and analytical methods can be used to determine the steady state resistance of the borehole. The numerical methods require much more computational effort but are generally more accurate than approximate analytical methods such as the Gu and O’Neal (a 1998) equivalent diameter method. The general equation for borehole resistance comes from summing the three resistances (fluid, pipe, and grout) between the fluid inside the Utube and the borehole wall as discussed in section 1.2.1. The fluid resistance is typically calculated with a convection correlation. The pipe resistance is determined as a cylindrical conductive resistance. The grout resistance is more difficult to determine, due to the complex geometry. 51 The methods that are compared in this chapter for calculating the grout resistance are the multipole method, the Paul (1996) method, the Gu and O’Neal (a 1998) approximate diameter method, and the cylinder source method (Ingersoll, 1948, 1954). The two numerical programs that are used to calculate the borehole resistance are GEMS2D (Rees, 2001) and Yavuzturk’s (Yavuzturk and Spitler, 1999) pie sector approximation which both use the finite volume method. 2.1 Borehole Resistance Transient and Steady State For the first few hours of constant heat injection or extraction the borehole resistance is transient. Figure 21 shows how the borehole resistance changes over the first 12 hours after a constant flux is applied. Figure 21 was generated from the average fluid temperature and the average temperature at the borehole wall radius using Equation 21. ( ) Q T T R f borehole total − = where, total R = borehole resistance ⎟⎠ ⎞ ⎜⎝ ⎛ W Km or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF f T = Average fluid temperature (K) or (ºF) borehole T = Average temperature at the borehole wall radius (K) or (ºF) Q = flux per unit length ⎟⎠ ⎞ ⎜⎝ ⎛ m W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft Btu (21) This figure comes from a GEMS2D simulation where the borehole was 11.4 cm (4.5 in) diameter, with a 1.6 cm (0.63 in) shank spacing, standard bentonite grout, a pipe 52 conductivity of 0.39 (W/(m·K)) (0.225 Btu/(h·ft·°F)), and a fluid convection coefficient of 1690 (W /(m2K)) (298 (Btu /(h⋅ ft2 ⋅°F))) As can be seen the borehole resistance is almost constant after about 12 hours. For typical boreholes there is usually less than a 2% difference between the steady state value and the value at 10 hours. This is indicated by Figure 22. Resistance vs Time 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0 1 2 3 4 5 6 7 8 9 10 11 12 Time (hours) Resistance (Km/W) Figure 21 Transient Borehole Resistance Profile vs Time Percent Difference of Transient Resistance with Respect to Steady State Resistance vs Time 0 10 20 30 40 50 60 70 80 90 100 0 1 2 3 4 5 6 7 8 9 10 11 12 Time (hours) Percent Difference (%) Figure 22 Percent Difference in Transient Borehole Resistance with respect to Steady State borehole resistance 53 The rate at which the resistance approaches the steady state value is dependent on the geometry and thermal properties within the borehole. A borehole with a low grout or pipe conductivity requires more time for the borehole to reach steady state. 2.2 Borehole Resistance Calculation from Analytical and Empirical Methods The analytical methods that are compared in this thesis include the multipole method, the Gu and O’Neal approximate diameter method, and the cylinder source method. The Paul method is also included in this section because it is based on curve fits to numerical and analytical data and is not strictly numerically based. The literature review in sections 1.2.1.1 through 1.2.1.4 describes how each method, except for the cylinder source, can be used to model the steady state borehole resistance. An application of the cylinder source for calculating the steady state borehole resistance is described in this section. The Gu and O’Neal method and the Paul method are much simpler to calculate than the multipole and cylinder source methods; however the multipole and cylinder source methods produce more accurate solutions. The Gu and O’Neal method was used exactly as described in section 1.2.1.2 and the Paul method was used exactly as described in section 1.2.1.3. Thus, with regards to these methods, no further explanation is necessary in this chapter. However, because of the complexity of the cylinder source and multipole methods, additional explanation is provided here in addition to what has been previously described in the literature review in sections 1.2.1.4 and 1.2.1.5 respectively. 54 The cylinder source solution can be used to model steady state resistance by using the principal of superposition. To model the grout resistance using Equation 14, the temperature rise from the flux exiting each leg of the Utube can be superimposed to calculate the average Utube outside wall temperature and also the average borehole wall temperature. Implicit in this method is the assumption that the soil conductivity, which is different from grout conductivity, has relatively little influence on the borehole thermal resistance. This method will be explained in greater detail using Figure 23 which shows the cylinder source locations for the calculation of the Utube average outside wall temperature. Figure 23 (a) shows a cylinder source located at (x,y) = (0,0) in an infinite medium of grout. The circle labeled borehole in Figure 23 (a) is shown to indicate the cylinder source location inside the borehole. The cylinder source in Figure 23 (a) is the location where the Utube temperature ( U tube T − ) will be calculated for Equation 14. A circle showing the borehole radius is shown, however the cylinder source method does not make a distinction between conductivities or thermal properties nor does it account for the existence of the other Utube since the medium is infinite. For the purpose of calculating borehole resistance, the conductivity of the soil ( soil k ) in Equation 113 should be replaced with the conductivity of the grout ( grout k ). In Figure 23 (a) the temperature rise at the U–tube radius should be calculated using r = o r in Equation 113. Since the temperature rise will not vary it is unnecessary to calculate several temperatures along the Utube circumference and average them. 55 Figure 23 Cylinder Source Diagram for Calculating the Utube Outside Wall Temperature for use in the Steady State Borehole Resistance Calculation Figure 23 (b) is similar to (a) except it models the cylinder source at the other leg of the Utube. Similar to Figure 23 (a), the cylinder source in (b) is also surrounded by an infinite medium of grout. To show where the temperature rise will be calculated, a circle is drawn showing Utube leg 1, however neither the fluid nor the Utube leg 1 pipe exists in the infinite and continuous medium surrounding the cylinder source. A larger 56 view of the Utubes in Figure 23 (b) is shown in Figure 23 (c) and the geometry of how to calculate the radius (r) for Equation 113 is shown in Figure 23 (d). The temperature rise from the cylinder source at Utube leg 2 should be calculated for several points along the circumference of Utube leg 1. The average temperature should then be calculated. This can be done by calculating “r” with the equation shown in Figure 23 (d) for several different Φ angles. To calculate the overall temperature rise at the Utube radius the average temperature increase from each Utube cylinder source should be superimposed to yield the overall temperature. Equation 25 calculates the resistance between the Utube OD and radius infinity. In a real system this is analogous to calculating the combined grout and soil resistance except the soil properties are the same as grout. Likewise Equation 26 calculates the resistance between borehole OD and radius infinity also using grout properties. Equation 27 finds the resistance of the grout that is located between the Utube and the borehole radius by subtracting the resistance calculated in Equation 26 from that calculated in Equation 25. In a similar manner the average temperature at the borehole radius can be found by averaging the temperature rises created by a cylinder source located at each Utube leg. If the borehole is symmetrical down the middle then the calculation shown in Equation 26 will calculate the resistance between the borehole wall and an infinite radius. 57 r (θ ) d 2 r 2 2r d cos(θ ) BH borehole borehole = + − r (θ ) (d / 2)2 r 2 2r (d / 2)cos(θ ) U tube UT UT = + − − ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ = ⋅ k r G F r k R r r k o o o ( , ,α , ) 1 , , Σ− = − − ⎥⎦ ⎤ ⎢⎣ = ⋅ ⎡ ⋅ + 1 0 ( ( 360), , , ) ( , , , ) 2 1 N U tube u tube UT grout grout UT UT grout grout r k R r r k N R r N R θ θ α α Σ− = ⎥⎦ ⎤ ⎢⎣ = ⋅ ⎡ ⋅ 1 0 1 ( ( 360), , , ) N BH BH UT grout grout r k N R r N R θ θ α grout U tube BH R = R − R − where, (22) (23) (24) (25) (26) (27) grout R = borehole grout resistance ⎟⎠ ⎞ ⎜⎝ ⎛ W mK or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF U tube R − = resistance between Utube OD and radius infinity ⎟⎠⎞ ⎜⎝ ⎛ W mK or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF BH R = resistance between borehole OD and radius infinity ⎟⎠ ⎞ ⎜⎝ ⎛ W mK or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF N = number of points to calculate the resistance around the circumference of the borehole or the Utube (dimensionless) θ = dummy variable used for counting from 0 to N (dimensionless) UT r = radius of the Utube (m) or (in) grout α = thermal diffusivity of the grout ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ s m2 or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ s ft 2 grout k = grout conductivity ⎟⎠ ⎞ ⎜⎝ ⎛ mk W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft⋅ F Btu o R(r, r , , k) o α = general function for resistance as a function of G() from Equation 115 (θ ) U tube r − = radial distance from the center of the cylinder source and the outside diameter of the other leg of the Utube (m) or (in) (θ ) BH r = radial distance from the center of the cylinder source and the outside diameter of the other leg of the Utube (m) or (in) d = the distance between the two centers of the left and right Utube legs (m) or (in) borehole r = borehole radius (m) or (in) 58 As shown in Equation 27, to calculate the grout resistance the resistance between the borehole wall and infinity is subtracted from the resistance between the Utube OD and infinity. Thus, Equations 22 through 27 can be applied to yield a solution for the steady state borehole resistance. In order to use the cylinder source to calculate steady state resistance a Fourier number should be calculated and the number of points to solve for around the borehole and Utube radiuses should be established. The time that was used to calculate the Fourier number was 80,000 hours which causes the solution to converge to five or more digits. The number of points that was chosen for finding the average temperature of both the Utube and borehole wall was 90. As the number of points was increased from 8 to180, the solution converged to five or more digits at 90 points. Other parameters that were chosen were the integration bounds in Equation 115 which are between zero and infinity. It was found that an upper integration bound of 10,000 produces 4 or more significant digits of convergence as compared to an integration bound of 100,000 which gives more than 8. The multipole resistance was found using a modified version of the Fortran 77 source code given in Bennet and Claesson (1987). Within the multipole method, the borehole resistance is found by establishing a temperature at the Utube wall and then calculating a heat flux and a temperature profile around the circumference of the borehole wall. The temperature at the borehole was calculated by taking an average of 180 points along the circumference of the borehole wall. Averaging 180 points versus averaging 360 points produced a temperature difference of less than 0.00001 °C (0.000018 °F) difference. The resistance can be calculated by using Equation 11. 59 2.3 Borehole Resistance Calculation using Numerical Methods As described in section 1.2.2.3, GEMS2D closely approximates the borehole geometry using a boundary fitted grid, whereas Yavuzturk approximates the borehole geometry using a pie shaped wedge in a parametric grid to represent the Utubes. In this chapter, GEMS2D is used as a standard for the borehole resistance calculation since it correlates very closely with another highly accurate method, the multipole method. This will be shown in section 2.5. GEMS2D has also proven to be a very accurate two dimensional finite volume program for other simulations. The disadvantage of GEMS2D is that it is approximately half as fast as Yavuzturk’s finite volume model program. In both GEMS2D and Yavuzturk’s pie sector approximation the average borehole wall temperature was subtracted from the given fluid temperature. This is shown in Equation 21. With constant flux and large times, typically greater than 10 hours, Equation 21 produces the steady state borehole resistance. A comparison of the steady state borehole resistances that GEMS2D and Yavuzturk’s pie sector approximation produce is shown in section 2.4. 2.4 Numerical Methods: Comparison between GEMS2D and the PieSector Approximation for Calculating Steady State Resistance Table 21 shows the baseline borehole system configuration that was used for the comparison. By varying individual parameters, this borehole configuration produced Table 22 which shows the steady state borehole resistance for both GEMS2D and the pie sector approximation. The conductivities that were varied in this comparison are the soil, grout and pipe conductivity. 60 The conductivities that were shown cover most typical borehole configurations. Also the borehole diameters that were chosen are 11.4 cm (4.5 in), 15.2 cm (6 in) and 19.1 cm (7.5 in) which also cover typical borehole configurations. As the borehole diameter was changed the shank spacing was held constant at 0.16 cm (0.067 in). The final parameter that was varied was the fluid flow rate at 0.000189 m3/s (3 gpm), 0.000379 m3/s (6 gpm), and 0.000568 m3/s (9 gpm). This also covers the range of most boreholes. Table 21 Borehole Properties (Base Case) Borehole System Table English Units SI Units Diameter 4.5 (in) 114.3 (mm) Shank Spacing 0.067 (in) 1.7 (mm) Utube OD 1.05 (in) 26.67 (mm) Utube ID 0.824 (in) 20.93 (mm) soil k 1.2 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft⋅ F Btu o 2.077 ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W o grout k 0.4 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft⋅ F Btu o 0.692 ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W o pipe k 0.8 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft⋅ F Btu o 1.38 ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W o fluid k 0.8 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft⋅ F Btu o 1.38 ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W o Fluid Volumetric Heat 62.4 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ft ⋅ F Btu 3 o 4.18 ⎟⎠ ⎞ ⎜⎝ ⎛ m ⋅K MJ 3 Flow Rate 3 (gpm) 0.000189 (m3 / s) 61 Table 22 Borehole Resistance Comparison between GEMS2D and the pie sector approximation Varied Input Borehole Resistance % Input s English SI Pie Sector GEMS2D Diff. (oF / h⋅ Btu) (K /W) (oF / h⋅Btu) (K /W) % soil k 0.8 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h⋅ ft ⋅°F Btu 1.38 ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W o 0.3615 0.685 0.3603 0.683 0.333 soil k 1.2 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft ⋅°F Btu 2.07 ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W o 0.3608 0.684 0.3588 0.680 0.559 soil k 1.6 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft ⋅°F Btu 2.77 ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W o 0.3604 0.683 0.3580 0.679 0.671 grout k 0.2 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft ⋅°F Btu 0.346 ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W o 0.6778 1.28 0.6481 1.23 4.48 grout k 0.4 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h⋅ ft ⋅°F Btu 0.692 ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W o 0.3608 0.684 0.3588 0.680 0.559 grout k 0.8 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h⋅ ft ⋅°F Btu 1.38 ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W o 0.1982 0.376 0.2143 0.406 7.79 pipe k 0.4 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h⋅ ft ⋅°F Btu .692 ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W o 0.3898 0.739 0.4284 0.812 9.43 pipe k 0.8 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h⋅ ft ⋅°F Btu 1.38 ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W o 0.3608 0.684 0.3588 0.680 0.559 pipe k 1.2 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h⋅ ft ⋅°F Btu 2.07 ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W o 0.3494 0.662 0.3329 0.631 4.84 Dia. 4.5 (in) 11.4 (cm) 0.3608 0.684 0.3588 0.680 0.559 Dia. 6 (in) 15.2 (cm) 0.4259 0.807 0.3138 0.595 30.3 Dia. 7.5 (in) 19.1 (cm) 0.4758 0.902 0.2785 0.528 52.3 Flow Rate 3 (gpm) 0.000189 (m3 / s) 0.3608 0.684 0.3588 0.680 0.559 Flow Rate 6 (gpm) 0.000379 (m3 / s) 0.3585 0.680 0.3570 0.677 0.413 Flow Rate 9 (gpm) 0.000568 (m3 / s) 0.3576 0.678 0.3563 0.675 0.358 Several observations can be made between the steady state resistance obtained from GEMS2D and the pie sector approximation in Table 22. The pie sector approximation deviates from the GEMS2D solution when the grout geometric properties are changed. By changing the diameter of the borehole without changing the shank spacing or the Utube diameter the grout geometry is being changed. This test measures 62 how accurately the automated grid generation algorithm approximates the actual geometry of the borehole with the pie sector approximation. The resistance is overestimated by more than 50% in the 19.1 cm (7.5 in) diameter case. When the conductivity of the pipe or grout are changed from the standard of 1.38 and 0.692 ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W o (0.8 and 0.4 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h⋅ ft ⋅°F Btu ) respectively the relative percent difference increases substantially. This error is accounted for because, when the conductive properties of the grout and pipe are increased or decreased, the effects of the geometric differences between the two programs are magnified. Changing the conductivity of the soil does not appreciably change the percent difference between the two programs. This is because soil conductivity has a second order effect on borehole resistance since it is outside the borehole and both programs accurately represent the circular geometry at the borehole wall radius. Both GEMS2D and the pie sector approximation would be poor choices as the resistance calculator for use with the BFTM model since they are very slow. The pie sector approximation requires on the order of thirty minutes and GEMS2D requires an hour on a 450 MHZ computer. Therefore analytical methods need to be compared to arrive at a reasonable solution. The pie sector approximation will not be considered further since GEMS2D is the more accurate finite volume model program for predicting borehole resistance as shown in section 2.5. 2.5 Comparison of Methods for Calculating Steady State Borehole Resistance The steady state borehole resistance calculation methods that are compared in this chapter include the Paul (1996) method, the Gu and O’Neal (a 1998) approximate 63 diameter method, cylinder source, and the multipole methods. The GEMS2D solution is given as a general comparison. The data for the baseline borehole used in this study are given in Table 23. Two different grout types were chosen, standard bentonite and thermally enhanced grout. This will give an understanding for how grout conductivity affects the different borehole resistance calculation methods. For each grout type, resistances were calculated for three different borehole diameters 7.6 cm (3 in), 11.4 cm (4.5 in), and 15.2 cm (6 in). The 7.6 cm (3 in) diameter case is an unrealistic borehole configuration; however it is useful for testing the capabilities of the models. For each borehole diameter, resistances were calculated for four different Utube shank spacings ranging from 3.2 mm (0.125 in) from the outside wall of each Utube, to where both Utubes are touching the borehole outside wall. The parameters that were not varied include the Utube diameter, Utube thermal properties, soil thermal properties, and the circulating fluid’s convection coefficient. Table 23 Base Line Borehole Properties Borehole Diameter Pipe  1" SDR11 D = 114.30 mm 4.5 in I.D. = 27.4 mm 1.08 in Grout  Standard Bentonite O.D. = 33.4 mm 1.31 in K = 0.75 ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W o .433 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft ⋅°F Btu K = 0.390 ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W o 0.225 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft ⋅°F Btu p C ρ = 3.90 ⎟⎠ ⎞ ⎜⎝ ⎛ m K MJ 3 58.2 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ft ⋅°F Btu 3 p C ρ = 1.77 ⎟⎠ ⎞ ⎜⎝ ⎛ m K MJ 3 26.4 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ft ⋅°F Btu 3 Soil  Typical Properties Spacing K = 2.50 ⎟⎠⎞ ⎜⎝ ⎛ m⋅ K W o 1.44 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h⋅ ft ⋅°F Btu A S1= 3.18 mm 1/8 in p C ρ = 2.50 ⎟⎠ ⎞ ⎜⎝⎛ m K MJ 3 37.3 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ft ⋅°F Btu 3 B S1= S2 S2 Fluid Convection Coefficient C3 S2= 3.00 mm 0.12 in H = 1690 ⎟⎠ ⎞ ⎜⎝ ⎛ m K W 2 298 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft ⋅°F Btu 2 C S2= 0 mm 0 in S1 S2 64 Table 24 gives the resistances that were calculated using the given borehole properties with the different methods. Table 25 shows the percent error of the steady state borehole resistance with respect to the GEMS2D calculated borehole resistance. Since GEMS2D is not capable of calculating the borehole resistance for the C spacing case, the multipole solution was used for the error calculation. The general methods for calculating the borehole resistance were shown in the literature review for the cylinder source, the multipole, Gu and O’Neal, and the Paul methods. The cylinder source column shown in Table 24 shows data using the cylinder source method described in section 2.2. The multipole data shown in Table 24 uses the tenth order multipole solution. 65 Table 24 Steady State Borehole Resistance Comparison N/A signifies that the method was not suitable for calculating the borehole resistance for this case Borehole Diameter Utube Spacing grout k Paul Gu and O'Neal Cylinder Source Multipole GEMS2D (mm) (in) (W / Km) ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅° Btu hr ft F (Km/W) ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅° Btu hr ft F (Km/W) ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅° Btu hr ft F (Km/W) ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅° Btu hr ft F (Km/W) ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅° Btu hr ft F (Km/W) ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅° Btu hr ft F 76.2 (3) A 0.75 (0.4333) 0.188 (0.326) 0.136 (0.235) 0.1337 (0.2314) 0.1213 (0.2099) 0.1211 (0.2095) 76.2 (3) B 0.75 (0.4333) 0.170 (0.294) 0.136 (0.235) 0.1338 (0.2316) 0.1214 (0.2101) 0.1213 (0.2099) 76.2 (3) C3 0.75 (0.4333) N/A 0.135 (0.233) 0.1331 (0.2303) 0.1206 (0.2087) 0.1204 (0.2084) 76.2 (3) C 0.75 (0.4333) 0.127 (0.220) 0.119 (0.206) 0.1170 (0.2025) 0.1025 (0.1774) N/A 114.3 (4.5) A 0.75 (0.4333) 0.256 (0.443) 0.222 (0.384) 0.2197 (0.3803) 0.2119 (0.3668) 0.2116 (0.3663) 114.3 (4.5) B 0.75 (0.4333) 0.205 (0.354) 0.190 (0.329) 0.1882 (0.3258) 0.1823 (0.3155) 0.1822 (0.3154) 114.3 (4.5) C3 0.75 (0.4333) N/A 0.146 (0.252) 0.1437 (0.2487) 0.1288 (0.2230) 0.1288 (0.2230) 114.3 (4.5) C 0.75 (0.4333) 0.141 (0.244) 0.137 (0.238) 0.1355 (0.2346) 0.1149 (0.1989) N/A 152.4 (6) A 0.75 (0.4333) 0.322 (0.557) 0.283 (0.489) 0.2807 (0.4859) 0.2737 (0.4737) 0.2734 (0.4733) 152.4 (6) B 0.75 (0.4333) 0.235 (0.407) 0.227 (0.392) 0.2248 (0.3892) 0.2216 (0.3836) 0.2216 (0.3836) 152.4 (6) C3 0.75 (0.4333) N/A 0.163 (0.282) 0.1611 (0.2788) 0.1386 (0.2399) 0.1387 (0.2401) 152.4 (6) C 0.75 (0.4333) 0.152 (0.263) 0.157 (0.273) 0.1556 (0.2694) 0.1260 (0.2182) N/A 76.2 (3) A 1.5 (0.8666) 0.116 (0.201) 0.0896 (0.155) 0.08779 (0.1520) 0.08796 (0.1523) 0.08774 (0.1519) 76.2 (3) B 1.5 (0.8666) 0.107 (0.185) 0.0896 (0.155) 0.08786 (0.1521) 0.08803 (0.1524) 0.08788 (0.1521) 76.2 (3) C3 1.5 (0.8666) N/A 0.0893 (0.155) 0.08743 (0.1513) 0.08763 (0.1517) 0.08740 (0.1513) 76.2 (3) C 1.5 (0.8666) 0.0853 (0.148) 0.0813 (0.141) 0.07947 (0.1376) 0.07899 (0.1367) N/A 114.3 (4.5) A 1.5 (0.8666) 0.150 (0.259) 0.133 (0.230) 0.1308 (0.2264) 0.1317 (0.2280) 0.1315 (0.2276) 114.3 (4.5) B 1.5 (0.8666) 0.124 (0.215) 0.117 (0.202) 0.1151 (0.1992) 0.1158 (0.2005) 0.1157 (0.2003) 114.3 (4.5) C3 1.5 (0.8666) N/A 0.0946 (0.164) 0.09279 (0.1606) 0.09149 (0.1584) 0.09144 (0.1583) 114.3 (4.5) C 1.5 (0.8666) 0.0922 (0.160) 0.0905 (0.157) 0.08871 (0.1536) 0.08627 (0.1493) N/A 152.4 (6) A 1.5 (0.8666) 0.183 (0.316) 0.163 (0.282) 0.1613 (0.2793) 0.1624 (0.2811) 0.1621 (0.2806) 152.4 (6) B 1.5 (0.8666) 0.139 (0.241) 0.135 (0.234) 0.1334 (0.2309) 0.1345 (0.2328) 0.1344 (0.2327) 152.4 (6) C3 1.5 (0.8666) N/A 0.103 (0.179) 0.1015 (0.1757) 0.09828 (0.1701) 0.09833 (0.1702) 152.4 (6) C 1.5 (0.8666) 0.0978 (0.169) 0.101 (0.174) 0.09876 (0.1710) 0.09413 (0.1629) N/A 66 Tables 24 and 5 show that the tenth order multipole and GEMS2D correlate very closely yielding a maximum difference of 0.26 percent. In addition, the multipole method is very fast with a computer compared to GEMS2D. It takes less than a second to calculate on a 450 MHz computer whereas the GEMS2D program might require half an hour, depending on the grid, on the same computer. Table 25 Percent Error of Borehole Resistance N/A signifies that the method was not suitable for calculating the borehole resistance for this case Bore Diameter Spacing Kgrout Paul Gu and O'Neal Cylinder Source Multipole mm (in) (W / Km) ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅° Btu hr ft F 76.2 (3) A 0.75 (0.4333) 55.5 11.9 10.4 0.17 76.2 (3) B 0.75 (0.4333) 39.9 11.8 10.4 0.11 76.2 (3) C3 0.75 (0.4333) N/A 12.0 10.5 0.17 76.2 (3) C 0.75 (0.4333) 23.8 15.9 14.1 N/A 114.3 (4.5) A 0.75 (0.4333) 20.8 4.69 3.83 0.12 114.3 (4.5) B 0.75 (0.4333) 12.3 4.28 3.29 0.03 114.3 (4.5) C3 0.75 (0.4333) N/A 12.9 11.5 0.01 114.3 (4.5) C 0.75 (0.4333) 22.5 19.5 17.9 N/A 152.4 (6) A 0.75 (0.4333) 17.8 3.35 2.67 0.09 152.4 (6) B 0.75 (0.4333) 6.14 2.31 1.46 0.01 152.4 (6) C3 0.75 (0.4333) N/A 17.5 16.1 0.10 152.4 (6) C 0.75 (0.4333) 20.6 24.9 23.5 N/A 76.2 (3) A 1.5 (0.8666) 32.2 2.10 0.05 0.25 76.2 (3) B 1.5 (0.8666) 21.3 2.00 0.02 0.17 76.2 (3) C3 1.5 (0.8666) N/A 2.13 0.04 0.26 76.2 (3) C 1.5 (0.8666) 7.96 2.86 0.60 N/A 114.3 (4.5) A 1.5 (0.8666) 13.9 0.86 0.50 0.20 114.3 (4.5) B 1.5 (0.8666) 7.25 0.95 0.59 0.07 114.3 (4.5) C3 1.5 (0.8666) N/A 3.43 1.48 0.05 114.3 (4.5) C 1.5 (0.8666) 6.89 4.89 2.82 N/A 152.4 (6) A 1.5 (0.8666) 12.8 0.64 0.48 0.17 152.4 (6) B 1.5 (0.8666) 3.72 0.57 0.77 0.02 152.4 (6) C3 1.5 (0.8666) N/A 5.03 3.22 0.05 152.4 (6) C 1.5 (0.8666) 3.91 6.81 4.92 N/A Compared to the tenth order multipole method, the Gu and O’Neal approximate diameter method, the cylinder source methods as well as the Paul method typically have greatly reduced accuracy. For all of the models in Table 25, the largest errors occur for 67 the 7.6 cm (3 in) diameter Utube case for both thermally enhanced and non thermally enhanced grout. The Paul method, the Gu and O’Neal method, and the cylinder source method all tend to over predict borehole resistance. The cylinder source is in some cases higher and some cases lower than the actual resistance as can be seen in Table 24. For both the Gu and O’Neal method and the cylinder source method, as the shank spacing increases from “A” (narrowly spaced Utube) to “C” (widely spaced Utube) the error increases substantially. For the Gu and O’Neal method, with standard grout, 11.4 cm (4.5 in) diameter with the “A” and “C” shank spacing, borehole errors were 4.7% and 19.5% respectively. For the same condition the cylinder source solution produced errors of 3.8% for the “A” shank spacing and 17.9% for the “C” shank spacing. This increase in error stems from the Gu and O’Neal method and the cylinder source method not taking into account the soil conductivity. As the Utubes move from very close together to very far apart the impact of soil conductivity on the borehole resistance increases. Thus, as would be expected for the thermally enhanced grout cases, the errors have all substantially decreased for both the Gu and O’Neal and the cylinder source methods, due to the grout and the soil conductivities being closer together. As stated earlier, the data in Table 25 shows that the Gu and ONeal method has an increase in error as the shank spacing increases. Thus, since the resistance is a direct result of the equivalent diameter (Equation 111), the data shows that the equivalent diameter calculation is less accurate for large shank spacings versus small shank spacings. The Paul method performed poorly in comparison to the other methods. In most cases the error produced by the Paul method was several times that of the other methods 68 shown in Table 25. Also, the error fluctuates differently with shank spacing than the Gu and O’Neal method and the cylinder source model. As mentioned in section 1.2.1.3 the experimental model had uniform heat flux around the Utubes which is not the case in a real system. This is the cause of some of the error in the Paul method however it probably does not account for all of the error in the 76.2 mm (3 in) diameter cases shown in Table 2.5. 2.6 Borehole Resistance and Merging of the Short and Long Time Step GFunction The steady state borehole resistance parameter is used to separate the long time step gfunction from specific borehole geometries making a single long time step gfunction valid for any specific borehole geometry. This is accomplished by the total R term in Equation 28. ( ) Q k T R Q T g soil f total ground − ⋅ − = 2π ( ) where, g = gfunction (nondimensionalized) total R = borehole resistance ⎟⎠ ⎞ ⎜⎝ ⎛ W Km or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF f T = average fluid temperature (K) or (ºF) ground T = steady state ground temperature (K) or (ºF) Q = flux per unit length ⎟⎠ ⎞ ⎜⎝ ⎛ m W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft Btu soil k = soil conductivity ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ft⋅ F Btu o (28) 69 The gfunction in Equation 28 is only a representation of the thermal resistance of the ground. Before the long time step gfunction can be used to calculate the fluid temperature the borehole resistance must be calculated using specific borehole parameters and then added to the thermal resistance of the ground. If the resistance calculation is not accurate then the long and short time step gfunctions will merge poorly. Figure 21 shows a short time step gfunction calculated with the line source for the borehole with properties shown in Table 21 with a 11.4 cm (4.5 in) diameter and B shank spacing. The steady state borehole resistance is shown in Table 24. In Figure 21 the long time step gfunction is for a single borehole. As can be seen in Figure 21 three different curves have been created for the longtime step gfunction using three different methods for calculating borehole resistance. The “LTS: Generalized” curve is the long time step gfunction without the borehole resistance. 70 Line Source Short Time Step Gfunction Compared to Long Time Step Gfunction Translated Using Different Borehole Resistance Calculation Methods 1 0 1 2 3 4 5 6 7 8 14 13 12 11 10 9 8 7 6 5 4 3 2 1 0 1 2 3 Log Time (ln(t/ts)) gfunction STS: Line Source LTS: GEMS2D and Multipole LTS: GuO'Neal LTS: Paul LTS: Generalized Figure 24 Line Source STS Gfunction Compared to LTS Gfunction Using Different Borehole Resistance Calculation Methods for a Single Borehole System As can be seen in Figure 23 the LTS and STS gfunction merges well using the resistance calculated with either the GEMS2D or Multipole resistance methods. Also the LTS gfunction using the Gu and O’Neal or the Paul methods matches less well with the STS gfunction. As shown in Table 25, the errors in the borehole resistances are 12.3% for the Paul method and 4.3% for the Gu and O’Neal method. The percent errors shown for this particular case in Table 25 are not the greatest errors. For some cases the merging between the long and short time step gfunctions will be even worse using the GuO’Neal and the Paul methods. 71 2.7 Conclusion As discussed in the literature review the long and short time step gfunctions are produced using different methods. The long time step gfunction is produced using superposition with data from a two dimensional radialaxial finite difference model. The short time step gfunction is produced with a two dimensional analytical or experimental model of the cross section of the borehole. Before the gfunction can be used in a simulation, consistency must be checked between the two methods that produce the short and long time step gfunctions and borehole resistance. If the short and long time step gfunction do not merge well together this is evidence of a problem with the borehole resistance calculation or with the short or long time step gfunction itself. This study shows that since the Paul method, for most geometries, does not accurately calculate the borehole resistance and therefore does not ensure a good merge of the long and short time step gfunction, it should not be used in simulations. The Gu O’Neal method is superior to the Paul method and might be suitable in a simulation when a very simple method is needed. The user should be aware of the errors involved with this simple calculation as shown in Table 25. Of the methods that are compared in this chapter, the multipole method is the best analytical method for the purpose of merging the long and short time step gfunction. Also, since the borehole resistance for most simulations will only be computed once, for a given simulation, it is not necessary for the resistance calculator to be exceptionally fast. However using the finite volume methods such as the pie sector approximation or GEMS2D which require fifteen minutes and 30 minutes, respectively, on a 1.4 Ghz computer is not practical. Since the multipole method requires less then a second to calculate on a 450 Mhz Pentium II and attains a 72 very good correlation with the GEMS2D model it is a very good choice for the borehole resistance calculator. 73 3 SHORT TIME STEP GFUNCTION CREATION AND THE BOREHOLE FLUID THERMAL MASS MODEL (BFTM) The short time step gfunction can be generated by any program or equation that is capable of approximating a transient borehole fluid temperature profile over time. The simplest and fastest method for use in a comp
Click tabs to swap between content that is broken into logical sections.
Rating  
Title  Development, Verification, and Design Analysis of the Borehole Fluid Thermal Mass Model for Approximating Short Term Borehole Thermal Response 
Date  20041201 
Author  Young, Thomas Ray 
Keywords  Borehole, GLHE, Ground Loop, BFTM, Borehole Fluid Thermal Mass, Borehole Fluid 
Department  Mechanical Engineering 
Document Type  
Full Text Type  Open Access 
Abstract  The borehole fluid thermal mass (BFTM) model, described in this thesis, is an analytical model for the heat conduction in and around a boreholetype ground loop heat exchanger (GLHE). The advantage of the BFTM model over other analytical models is its ability to model the thermal mass of the circulating fluid. The thesis also applies the BFTM model to generating a nondimensional representation of the shorttime soil thermal resistance, called a gfunction, which is used in the GLHE design software, GLHEPRO. The changes in required depth due to changing internal borehole properties are determined for peakload and nonpeakload dominant buildings by performing 1440 sizing simulations in GLHEPRO. The performance of the BFTM model on an hourly basis was determined by comparing the BFTM model to a line source model using 6 hourly HVACSIM+ simulations. Six different methods for calculating the steady state borehole resistance are also reviewed in this thesis. 
Note  Thesis 
Rights  © Oklahoma Agricultural and Mechanical Board of Regents 
Transcript  DEVELOPMENT, VERIFICATION, AND DESIGN ANALYSIS OF THE BOREHOLE FLUID THERMAL MASS MODEL FOR APPROXIMATING SHORT TERM BOREHOLE THERMAL RESPONSE by THOMAS RAY YOUNG Bachelor of Science Oklahoma State University Stillwater, Oklahoma 2001 Submitted to the Faculty of the Graduate College of the Oklahoma State University in partial fulfillment of the requirements for the degree of MASTERS OF SCIENCE Oklahoma State University December, 2004 ii DEVELOPMENT, VERIFICATION, AND DESIGN ANALYSIS OF THE BOREHOLE FLUID THERMAL MASS MODEL FOR APPROXIMATING SHORT TERM BOREHOLE THERMAL RESPONSE Thesis Approved: Dr. Spitler Thesis Advisor Dr. Delahoussaye Dr. Fisher Dr. Emsile Dean of Graduate College iii ACKNOWLEDGEMENTS In order to give credit where it is due, my conscience requires that, at the very least, I mention God who sent his son Jesus Christ to die in my place, for my sins, and gives eternal life to everyone who believes in Jesus. Without the grace, mercy, strength and wisdom God has provided me in various ways, this thesis would never have been finished. He deserves ALL the honor and glory. Also I am indebted to Dr. Spitler for his willingness to take me on as one of his students. His expertise and experience were invaluable and his willingness to continue working with me after I took a job and left Stillwater is very much appreciated I would like to thank Dr. Delahoussaye for his mentoring early on in my research. Dr. Delahoussaye was essential in helping me gain the programming skills and practical understanding that I needed to succeed. To Aditya, Hayder, and Liu I appreciate your friendships, learned from your expertise, and enjoyed working with you greatly. I would especially like to thank Xiaowei for the hours he spent running simulation and examining my work. I would also like to thank my wife Rachel who provided moral support through a listening ear and plenty of cookies and brownies. I would like to thank my parents for instilling in me moral values and a good work ethic. Also, thank you for your prayers and financial support. iv I should also mention my Sunday school class for all there prayers and moral support. Hardly a Sunday went by without someone asking me how the thesis was going. v TABLE OF CONTENTS 1 INTRODUCTION ....................................................................................................1 1.1 BACKGROUND....................................................................................................................... 5 1.2 LITERATURE REVIEW......................................................................................................... 5 1.2.1 Steady State Modeling of Boreholes ....................................................................................... 7 1.2.1.1. Line Source Model..................................................................................................................... 11 1.2.1.2. GuO’Neal Equivalent Diameter Model .................................................................................... 15 1.2.1.3. Paul Model................................................................................................................................ 17 1.2.1.4. Cylinder Source Model .............................................................................................................. 19 1.2.1.5. Multipole................................................................................................................................... 21 1.2.2 Transient Modeling of Borehole Heat Exchangers............................................................... 31 1.2.2.1. Buried Electrical Cable Model................................................................................................... 32 1.2.2.2. Gfunction Model: Long Time Step........................................................................................... 35 1.2.2.3. Gfunction Model: Short Time Step........................................................................................... 41 1.2.3 GLHEPRO Version 3 Design Tool ....................................................................................... 47 1.3 OBJECTIVES......................................................................................................................... 48 2 COMPARISON OF BOREHOLE RESISTANCE CALCULATION METHODS...................................................................................................................... 50 2.1 BOREHOLE RESISTANCE TRANSIENT AND STEADY STATE ..................................................... 51 2.2 BOREHOLE RESISTANCE CALCULATION FROM ANALYTICAL AND EMPIRICAL METHODS... 53 2.3 BOREHOLE RESISTANCE CALCULATION USING NUMERICAL METHODS ............................... 59 2.4 NUMERICAL METHODS: COMPARISON BETWEEN GEMS2D AND THE PIESECTOR APPROXIMATION FOR CALCULATING STEADY STATE RESISTANCE..................................................... 59 2.5 COMPARISON OF METHODS FOR CALCULATING STEADY STATE BOREHOLE RESISTANCE.. 62 2.6 BOREHOLE RESISTANCE AND MERGING OF THE SHORT AND LONG TIME STEP GFUNCTION 68 2.7 CONCLUSION............................................................................................................................ 71 3 SHORT TIME STEP GFUNCTION CREATION AND THE BOREHOLE FLUID THERMAL MASS MODEL (BFTM)............................................................. 73 3.1 BOREHOLE FLUID THERMAL MASS MODEL ........................................................................... 74 3.2 GROUT ALLOCATION FACTOR USED TO IMPROVE ACCURACY ............................................. 77 3.3 FLUID MULTIPLICATION FACTOR IN THE BFTM MODEL....................................................... 78 3.4 IMPLEMENTATION OF THE BFTM MODEL.............................................................................. 79 3.4.1 Bessel Function Evaluation .................................................................................................. 80 3.4.2 BFTM Model  Solving the Integral...................................................................................... 81 3.4.3 Incorporating the Fluid Thermal Mass Model in a Design Program................................... 83 3.5 IMPROVING THE BFTM MODEL FOR SMALL TIMES USING LOGARITHMIC EXTRAPOLATION 86 3.5.1 Implementing Logarithmic Extrapolation............................................................................. 87 4 NUMERICAL VALIDATION OF THE BOREHOLE FLUID THERMAL MASS MODEL USING GEMS2D................................................................................ 90 vi 4.1 GEMS2D SIMULATIONS .......................................................................................................... 93 4.2 FINDING THE GROUT ALLOCATION FACTOR .......................................................................... 95 4.3 BOREHOLE DIAMETER VALIDATION WITH LINE SOURCE COMPARISON ............................ 100 4.4 SHANK SPACING VALIDATION ............................................................................................... 105 4.5 GROUT CONDUCTIVITY VALIDATION.................................................................................... 110 4.6 SOIL CONDUCTIVITY VALIDATION ........................................................................................ 112 4.7 GROUT VOLUMETRIC HEAT CAPACITY VALIDATION .......................................................... 116 4.8 BFTM MODEL FLUID FACTOR VALIDATION WITH GEMS2D............................................. 120 4.9 IMPLEMENTATION AND VALIDATION OF THE BFTME MODEL.......................................... 124 4.10 CONCLUSION OF BFTM MODEL VALIDATION...................................................................... 129 5 THE EFFECT OF THE BFTM MODEL ON GLHE DESIGN....................... 130 5.1 TEST BUILDINGS..................................................................................................................... 130 5.1.1 Church ............................................................................................................................... 130 5.1.2 Small Office Building.......................................................................................................... 132 5.1.3 Annual Loading .................................................................................................................. 133 5.2 GLHE DESIGN PROCEDURES................................................................................................. 141 5.3 SIMULATION RESULTS............................................................................................................ 147 5.3.1 Fluid Factor Results ........................................................................................................... 148 5.3.2 UTube Shank Spacing Results........................................................................................... 159 5.3.3 Borehole Diameter Results ................................................................................................. 164 5.1.4 Discussion of GLHEPRO Sizing Results ............................................................................ 168 6 HOURLY SIMULATION USING THE BFTM MODEL................................ 170 6.1 HVACSIM+ HOURLY SIMULATION...................................................................................... 171 6.2 LINE SOURCE AND BFTM MODEL COMPARISON USING A DETAILED HVACSIM+ MODEL 173 6.3 INFLUENCE OF THE FLUID MULTIPLICATION FACTOR ON SYSTEM DESIGN........................ 181 6.3.1 Fluid Factor Analysis with HVACSIM+ Simulation Tools................................................. 182 6.4 HVACSIM+ AND GLHEPRO COMPARISON........................................................................ 185 6.5 CONCLUSION.......................................................................................................................... 193 7 CONCLUSIONS AND RECOMMENDATIONS.............................................. 194 7.1 CONCLUSIONS........................................................................................................................ 194 7.2 RECOMMENDATIONS .............................................................................................................. 196 vii LIST OF TABLES Table Page Table 11 Paul Curve Fit Parameters used to Calculate the Steady State Grout Resistance .................................................................................................................................. 18 Table 12 Variable Input List for the Multipole Method.................................................. 23 Table 21 Borehole Properties (Base Case) ...................................................................... 60 Table 22 Borehole Resistance Comparison between GEMS2D and the pie sector approximation ........................................................................................................... 61 Table 23 Base Line Borehole Properties ......................................................................... 63 Table 24 Steady State Borehole Resistance Comparison................................................ 65 Table 25 Percent Error of Borehole Resistance............................................................... 66 Table 31 Borehole Properties Table ................................................................................ 75 Table 41 Borehole Properties for GEMS2D to BFTM Model Comparisons. ................. 91 Table 42 Borehole Diameter vs Time for Slope Matching. ............................................ 98 Table 43 GAF Dependent on Borehole Diameter, Shank Spacing and Fluid Factor...... 99 Table 44 Borehole Resistances and GAF for Diameter Validation Tests ..................... 101 Table 45 Borehole Resistances for Shank Spacing Validation Tests............................ 106 Table 46 Borehole Resistances for Shank Spacing Validation Tests............................ 110 Table 47 Borehole Resistances for Soil Conductivity Validation Tests ....................... 113 Table 48 NonDimensional Borehole Diameter vs Time for Slope Matching.............. 126 Table 49 GAF Dependent on NonDimensional Borehole Diameter, NonDimensional Shank Spacing and Fluid Factor ............................................................................. 126 Table 51 Church Building Description.......................................................................... 131 Table 52 Church Building Load Table for Different Locations.................................... 134 Table 53 Building Load Table for the Small Office Building....................................... 138 Table 54 GLHE Properties for the Church and Small Office Building......................... 141 Table 55 Undisturbed Ground Temperature Table for Various Cities.......................... 142 Table 56 Percent Change in GLHEPRO Sizing Depth with Respect to Varying Fluid Factor for PeakLoadDominant and NonPeakLoadDominant Buildings with Standard Grout ........................................................................................................ 154 Table 57 Percent Change in GLHEPRO Sizing Depth with Respect to Varying Fluid Factor for PeakLoadDominant and NonPeakLoadDominant Buildings with Thermally Enhanced Grout..................................................................................... 155 Table 58 Required Depth (ft) for Thermally Enhanced Grout, Calculated with GLHEPRO.............................................................................................................. 156 Table 59 Required Depth (ft) for Standard Grout, Calculated with GLHEPRO........... 157 viii Table 510 Percent Change in GLHEPRO Sizing Depth for Changes in Borehole Shank Spacing.................................................................................................................... 163 Table 511 Percent Change in GLHEPRO Sizing Depth for Changes in Borehole Diameter.................................................................................................................. 167 Table 61 Coefficients for the VS200 Climate Master Heat Pump ................................ 172 ix LIST OF FIGURES Figure Page Figure 11 Borehole system ................................................................................................ 2 Figure 12 Borehole Cross Section..................................................................................... 3 Figure 13 Crosssection of the Borehole and the Corresponding Thermal ΔCircuit (Hellström, 1991 p.78)................................................................................................ 7 Figure 14 Cross Section of a Borehole with Symmetry Line and the Corresponding Thermal Circuit........................................................................................................... 8 Figure 15 Actual Geometry vs Equivalent Diameter Approximation............................. 16 Figure 16 Types of shank spacing used in the Paul borehole resistance approximation. 18 Figure 17 Example of a 2D System for the Multipole Method. ...................................... 22 Figure 18 Source and Sink Location for a Single Pipe.................................................... 26 Figure 19 Steady State Temperature Field for a Borehole Heat Exchanger ................... 29 Figure 110 Temperature Change Around the Borehole Circumference.......................... 30 Figure 111 Diagram of a buried electrical cable and circuit ........................................... 33 Figure 112 Short and long time step gfunction without borehole resistance ................. 37 Figure 113 Short and long time step gfunctions with borehole resistance..................... 37 Figure 114 Twodimensional radialaxial mesh for a heat extraction borehole in the ground (Eskilson, 1987)............................................................................................ 39 Figure 115 Long Time Step gfunction for a 64 Borehole System in an 8x8 Configuration with Varying Borehole Spacing ........................................................ 41 Figure 116 Grid for a cross section of a borehole ........................................................... 43 Figure 117 Grid for PieSector Approximation (Yavuzturk, 1999) ................................ 45 Figure 21 Transient Borehole Resistance Profile vs Time .............................................. 52 Figure 22 Percent Difference in Transient Borehole Resistance with respect to Steady State borehole resistance........................................................................................... 52 Figure 23 Cylinder Source Diagram for Calculating the Utube Outside Wall Temperature for use in the Steady State Borehole Resistance Calculation .............. 55 Figure 24 Line Source STS Gfunction Compared to LTS Gfunction Using Different Borehole Resistance Calculation Methods for a Single Borehole System ............... 70 Figure 31 Average Fluid Temperature using the line source and GEMS2D model with fluid mass for a heat rejection pulse ......................................................................... 73 Figure 32 Integrated Function in the BEC Model ........................................................... 83 Figure 33 STS GFunction Translation ........................................................................... 85 Figure 34 BFTME, BFTM, and GEMS2D Fluid Temperature and Gfunction ............ 89 Figure 41 Borehole Geometries Simulated with GEMS2D ............................................ 92 Figure 42 Blocks for a Borehole System without Interior Cells for GEMS2D............... 93 x Figure 43 GEMS2D Grid of a Borehole with Soil .......................................................... 94 Figure 44 GEMS2D Grid of Utube and Fluid with 8 Annular Regions ........................ 95 Figure 45 Fluid Temperature vs Time of GEMS2D and BFTM with Varying GAF...... 96 Figure 46 Gfunction Slope vs GAF for the BFTM at 5 Hours ...................................... 97 Figure 47 Fluid Temperature and Gfunction for 15.24 cm (6 in) Diameter Borehole Using a Slope Matching Time of 3 Hours for (a) and (b) and 5 hours for (c) and (d) .................................................................................................................................. 98 Figure 48 Validation of the BFTM Model Using a GEMS2D Simulation with a 7.62 cm (3 in) Borehole ........................................................................................................ 102 Figure 49 Validation of the BFTM Model Using a GEMS2D Simulation with a 11.4 cm (4.5 in) Borehole ..................................................................................................... 103 Figure 410 Validation of the BFTM Model Using a GEMS2D Simulation with 15.2 cm (6 in) Borehole with 3.16 cm (1.24 in) Shank Spacing .......................................... 104 Figure 411 Validation of the BFTM Model Using a GEMS2D Simulation with a 19.1 cm (7.5 in) Borehole with a 4.12 cm (1.62 in) Shank Spacing..................................... 105 Figure 412 Validation of the BFTM Model Using a GEMS2D Simulation with a 0.316 cm (0.125 in) Shank Spacing ...................................................................................107 Figure 413 Validation of the BFTM Model Using a GEMS2D Simulation with a 2.25 cm (0.89 in) Shank Spacing .................................................................................... 108 Figure 414 Validation of the BFTM Model Using a GEMS2D Simulation with a 1 4.13 cm (5/8 in) Shank Spacing...................................................................................... 109 Figure 415 Validation of the BFTM Model Using a GEMS2D Simulation with Grout Conductivity of 0.25 W/(m·K) (0.144 Btu/(h·ft·°F)) .............................................. 111 Figure 416 Validation of the BFTM Model Using a GEMS2D Simulation with grout conductivity of 1.5 W/m·K (0.867 Btu/(h·ft·°F)).................................................... 112 Figure 417 Validation of the BFTM Model Using a GEMS2D Simulation With Soil Conductivity of 0.5 W/m·K (0.289 Btu/(h·ft·°F))................................................... 114 Figure 418 Validation of the BFTM Model Using a GEMS2D Simulation with a Soil Conductivity of 1.5 W/m·K (0.867 Btu/(h·ft·°F))................................................... 115 Figure 419 Validation of the BFTM Model Using a GEMS2D Simulation with a Soil Conductivity of 8 W/m·K (4.62 Btu/(h·ft·°F))........................................................ 116 Figure 420 Validation of the BFTM Model Using a GEMS2D Simulation with a Grout Volumetric Heat Capacity of 2 MJ/m3·k (29.8 Btu/ft3·F) ....................................... 118 Figure 421 Validation of the BFTM Model Using a GEMS2D Simulation with a Grout Volumetric Heat Capacity of 8 MJ/m3K (119 Btu/ft3·F)....................................... 119 Figure 422 Validation of the BFTM Model Using a GEMS2D Simulation with a 11.4 cm (4.5 in) BH Diameter 3 cm (1.18 in) Shank Spacing and 2 Times the Fluid.......... 121 Figure 423 Validation of the BFTM Model Using a GEMS2D Simulation with a 11.4 cm (4.5 in) BH Diameter 3 cm (1.18 in) Shank Spacing and 4 Times the Fluid.......... 122 Figure 424 Validation of the BFTM Model Using a GEMS2D Simulation with a 19.05 cm (7.5 in) BH Diameter 2.25 cm (0.886 in) Shank Spacing and 2 Times the Fluid ................................................................................................................................ 123 Figure 425 Validation of the BFTM Model Using a GEMS2D Simulation with a 19.05 cm (7.5 in) BH Diameter 2.25 cm (0.886 in) Shank Spacing and 4 Times the Fluid ................................................................................................................................ 124 xi Figure 426 Temperature Profile for the BFTME and GEMS2D Models with a 17.8 cm (7 in) Borehole Diameter and 4 cm (1.57 in) Shank Spacing Using an Interpolated GAF Value.............................................................................................................. 128 Figure 51 Monthly Church Heating Loads.................................................................... 135 Figure 52 Monthly Church Cooling Loads.................................................................... 136 Figure 53 Monthly Church Peak Heating Loads ........................................................... 137 Figure 54 Monthly Church Peak Cooling Loads........................................................... 137 Figure 55 Monthly Heating Loads for the Small Office Building ................................ 139 Figure 56 Monthly Cooling Loads for the Small Office Building ................................ 139 Figure 57 Monthly Peak Heating Loads for the Small Office Building........................ 140 Figure 58 Monthly Peak Cooling Loads for the Small Office Building ....................... 140 Figure 59 Raw Church Loads Single 2 Hour Peak Heat Load...................................... 144 Figure 510 Raw Church Loads for Birmingham AL..................................................... 144 Figure 511 One Peak of Hourly Loads for Tulsa Small Office Building...................... 145 Figure 512 One Work Week of Hourly Loads for Tulsa Small Office Building .......... 146 Figure 513 Short Time Step GFunction Comparison between the Line Source and the BFTM model with 0.1, 1, 2, 3 and 4 x Fluid Factor ............................................... 147 Figure 514 GLHEPRO Sized Depth vs Fluid Factor for a Church Building with a Borehole Diameter of 4.5 in (11.4 cm), Enhanced Grout, and Shank Spacing of 0.125 in (0.318 cm)................................................................................................. 149 Figure 515 GLHEPRO Sized Depth vs Fluid Factor for Church Building with a Borehole Diameter of 4.5 in (11.4 cm), Enhanced Grout, and Shank Spacing of 1.87 in (4.75 cm)............................................................................................................. 149 Figure 516 GLHEPRO Sized Depth vs Fluid Factor for a Church Building with a Borehole Diameter of 4.5 in (11.4 cm), Standard Bentonite Grout, and Shank Spacing of 0.125 in (0.318 cm)............................................................................... 150 Figure 517 GLHEPRO Sized Depth vs Fluid Factor for Church Building with a Borehole Diameter of 4.5 in (11.4 cm), Standard Bentonite Grout, and Shank Spacing of 1.87 in (4.75 cm)................................................................................... 151 Figure 518 GLHEPRO Sized Depth vs Fluid Factor for Small Office Building with a Borehole Diameter of 4.5 in (11.4 cm), Standard Bentonite Grout, and Shank Spacing of 0.125 in (0.318 cm)............................................................................... 152 Figure 519 GLHEPRO Sized Depth vs Fluid Factor for Small Office Building with a Borehole Diameter of 4.5 in (11.4 cm), Thermally Enhanced Grout, and Shank Spacing of 1.87 in (4.75 cm)................................................................................... 152 Figure 520 Borehole Resistance vs Shank Spacing for Church Building Using Standard Bentonite Grout and a Diameter of 4.5 in (11.4 cm) .............................................. 159 Figure 521 Sized Borehole Depth vs Shank Spacing for Church Building Using Standard Bentonite Grout, Diameter of 4.5 in (11.4 cm), Fluid Factor of 1.......................... 160 Figure 522 Sized Borehole Depth vs Shank Spacing for Church Building Using Standard Bentonite Grout, Diameter of 4.5 in (11.4 cm), Fluid Factor of 2.......................... 161 Figure 523 Sized Borehole Depth vs Shank Spacing for Office Building Using Standard Bentonite Grout, Diameter of 4.5 in (11.4 cm), Fluid Factor of 1.......................... 161 Figure 524 Sized Borehole Depth vs Shank Spacing for Office Building Using Standard Bentonite Grout, Diameter of 4.5 in (11.4 cm), Fluid Factor of 2.......................... 162 xii Figure 525 Borehole Resistance vs Diameter Using Thermally Enhanced Grout, Shank Spacing of 0.125 in (0.318 cm)............................................................................... 164 Figure 526 Sized Borehole Depth vs Diameter for Church Building Using Thermally Enhanced Grout, Shank Spacing of 0.125 in (0.318 cm), and Fluid Factor of 1.... 165 Figure 527 Sized Borehole Depth vs Diameter for Small Office Building Using Standard Bentonite Grout , Shank Spacing of 0.125 in (0.318 cm), and Fluid Factor of 1 ... 165 Figure 61 Three Component Model of a GLHE System for Hourly Loads .................. 171 Figure 62 Gfunction’s for Various Fluid Factors......................................................... 173 Figure 63 Detailed HVACSIM+ Model with Tulsa Loads for PeakLoadDominant Times....................................................................................................................... 175 Figure 64 Detailed HVACSIM+ Model with Tulsa Loads for PeakLoadDominant Times....................................................................................................................... 175 Figure 65 Detailed HVACSIM+ with Tulsa Loads....................................................... 176 Figure 66 Detailed HVACSIM+ with Tulsa Loads for a Short Duration Heat Pulse on a Long Duration Heat Pulse....................................................................................... 177 Figure 67 Detailed HVACSIM+ with Tulsa Loads for Long Duration Heat Pulses..... 178 Figure 68 Difference in Temperature between the LS and BFTM Models................... 179 Figure 69 Heat Pump Power Curve for the LS and BFTM Models.............................. 180 Figure 610 Heat Pump Power Curve for the LS and BFTM Models............................ 180 Figure 611 Detailed HVACSIM+ Model with Tulsa Loads ......................................... 181 Figure 612 GLHE Inlet Temperature for Different Fluid Factors................................. 183 Figure 613 GLHE Inlet Temperature for Different Fluid Factor .................................. 184 Figure 614 GLHE Inlet Temperature for Different Fluid Factors................................. 185 Figure 615 Typical Peak Loads for Small Office Building in Houston ........................ 186 Figure 616 Typical Peak Loads for Church Building in Nashville ............................... 187 Figure 617 Two Component GLHE Model................................................................... 187 Figure 618 Entering Temperature to the Heat Pump for a Church Building Located in Nashville, 1 x fluid, and 2 hour GLHEPRO Peak Duration ................................... 188 Figure 619 Entering Temperature to the Heat Pump for a Church Building Located in Nashville, 2 x fluid, and 2 hour GLHEPRO Peak Duration ................................... 189 Figure 620 Entering Temperature to the Heat Pump for a Church Building Located in Nashville, 4 x fluid, and 2 hour GLHEPRO Peak Duration ................................... 189 Figure 621 GLHEPRO and HVACSIM+ Maximum Yearly Entering Temperature to the Heat Pump for a Church Building Located in Nashville ........................................ 190 Figure 622 GLHEPRO and HVACSIM+ Minimum Yearly Entering Temperature to the Heat Pump for a Church Building Located in Nashville ........................................ 191 Figure 623 Entering Temperature to the Heat Pump for a Small Office Building Located in Houston, 1 x fluid, and 8 hour GLHEPRO Peak Duration................................. 192 Figure 624 Entering Temperature to the Heat Pump for a Small Office Building Located in Houston, 2 x fluid, and 8 hour GLHEPRO Peak Duration................................. 192 1 1 INTRODUCTION “As long as the earth endures, seedtime and harvest, cold and heat, summer and winter, day and night will never cease.” (Genesis 8:22, NIV) The basic needs of man which include keeping warm in winter and cool in summer have remained constant throughout history. The technology used to meet the needs has changed. People and animals have historically used caves and manmade holes as shelter from the elements. In this way humans have been extracting heat from the earth to keep warm in winter and using the earth to keep cool in summer for centuries. Modern man uses more refined methods for extracting and rejecting heat from the ground such as ground source heat pump (GSHP) systems. The term “ground source heat pump (GSHP)” refers to heat pump systems that use either the ground or a water reservoir as a heat source or sink. GSHP systems are either openloop or closedloop. Openloop GSHP systems use a pump to circulate groundwater through the heat pump heat exchanger. A closedloop GSHP system uses a water pump to circulate fluid through pipes buried horizontally or inserted into boreholes in the ground. The buried closed loop version of the GSHP is commonly referred to as a ground loop heat exchanger (GLHE). The physical properties of boreholes are very important to the study of GLHE systems. Boreholes typically range between 46 to 122 meters (150 to 400 ft) deep and are typically around 10 to 15 cm (4 to 6 inches) in diameter. A borehole system can be 2 composed of anywhere from 1 to over 100 boreholes. Each borehole in a multiborehole system is typically placed at least 4.5 m (15 ft) from all other boreholes. Figure 11 shows a vertical cross section of three boreholes. Each borehole is connected to the other boreholes with pipes that are typically buried 12 meters (36 ft) under the top surface. Figure 11 Borehole system After the Utube is inserted, the borehole will usually be backfilled with grout. The grout is used to prevent contamination of aquifers. Figure 11 shows 3 ideal boreholes. The grout is often a bentonite clay mixture, with the possibility of having thermallyenhanced additives. The grout usually has a thermal conductivity significantly lower than the surrounding ground. The circulating fluid is water or a waterantifreeze mixture. Each borehole is shown in this picture to be parallel with the other boreholes and perpendicular to the surface of the ground. In reality, drilling rigs do not drill perfectly straight, causing the path of a borehole to deviate, especially in deep boreholes. The Utube as shown in Figure 11 has equal spacing between the two legs of the Utube throughout the borehole. In real systems, however, the Utube leg spacing does 3 not necessarily remain constant throughout the length of the borehole. Spacers are sometimes employed to force the tubes towards the borehole wall. Figure 12 shows a two dimensional horizontal cross section of a single borehole. The Utube leg spacing is called the shank spacing and is defined as the shortest distance between the outer pipe walls of each leg of the Utube. UTUBE SHANK SPACING GROUT BOREHOLE GROUND Figure 12 Borehole Cross Section As previously mentioned, the size of a borehole heat exchanger system can range from one borehole to over a hundred. For small buildings one borehole may suffice but for large commercial buildings over 100 boreholes are sometimes required. This can make the initial investment quite costly. The main advantage of a GSHP system over an airsource heat pump system is that it rejects heat to the ground in the summer, when the ground is cooler than the air, and extracts heat from the ground during the winter, when the ground is warmer than the air. A GSHP system will very seldom reject the same amount of heat as it extracts on an annual basis. In cold climates, for envelopedominated buildings, the GSHP system will extract much more heat from the ground than it rejects to the ground. In this case, the ground surrounding the boreholes gradually declines in temperature. Over time, the reduction in the ground temperature around the boreholes will decrease the performance 4 of the heat pump in heating mode. In cold climates, the fluid circulating in the boreholes might drop below freezing, requiring the addition of antifreeze in the system. Similarly in warm climates, since more heat is rejected to the ground than extracted from the ground, the ground temperature will rise. This will impair the performance of the heat pump in cooling mode. The actual annual imbalance depends not only on climate but also on the building internal heat gains and building design. The thermal loads over a number of years must be accounted for when designing a GSHP system. This is necessary to determine the impact of any annual heat imbalance. If a borehole heat exchanger (BHE) is overdesigned, the initial construction costs may be excessive. If the system is underdesigned, the BHE may not meet the long term heating or cooling needs of the user. Research has been conducted for the purpose of applying GSHP technology to other areas besides buildings. Chiasson and Spitler (2001 and 2000) at Oklahoma State University have conducted research applying GSHP technology to highway bridges. The system uses pipes embedded in the road pavement to circulate fluid from a GSHP to eliminate ice or snow formation. The potential benefits of this new application include safer driving conditions and longer lasting bridges and roads due to reduced corrosion. Engineers who are attempting to design a GSHP system for a specific application can use programs such as GLHEPRO (Spitler 2000) to describe a potential system composed of a specific ground loop heat exchanger and heat pump and then simulate the systems response to monthly and peak, heating and cooling loads. Using programs such as GLHEPRO, engineers can also optimize the depth of a specific borehole heat exchanger configuration. 5 1.1 BACKGROUND Regardless of the GSHP application, the thermal response of the GLHE plays an important part in the design and simulation of GSHP systems. Since thermal loading on typical GLHE systems is of long duration, design methodologies have focused, in some detail, on long time step responses to monthly loads. (Eskilson, 1987) However, short time responses have typically been modeled only crudely using analytical models such as the cylinder source. These short time responses can be very important in determining the effect of daily peak loads. Daily peak loads occur in all applications but may be dominant in applications such as church buildings, concert halls, and the Smart Bridge application. To model the shorttime response, it is important to accurately represent such details as the borehole radius, Utube diameter and shank spacing, as well as the thermal properties and mass of the circulating fluid, Utube and grout. This thesis presents a new methodology for modeling the short time GLHE thermal response. This is particularly important for systems with peakloaddominant loading conditions. 1.2 LITERATURE REVIEW The literature review describes different methods that have been developed to model borehole heat exchangers. The methods are divided into two categories: steady state and transient. Quasisteady state conditions occur in twodimensional borehole cross sections, as shown in Figure 12 when the circulating fluid, Utubes and grout within a borehole do not change temperature (relative to each other) with time for a constant heat flux. If the internal borehole temperature differences are constant, the borehole resistance, defined as the resistance between the circulating fluid and the borehole, is also constant. Thus, 6 when a borehole’s internal temperature differences have stabilized for a constant heat flux, the borehole resistance can be modeled as a constant. Transient modeling of borehole heat exchangers might be broken into three different regions. The first region deals with transients that occur within the borehole before the borehole reaches steady state. For this transient region, the borehole may be modeled as having infinite length since surface and bottom end effects can be neglected. Two dimensional geometric and thermal properties of the borehole influence the temperature response in the first region. The second region occurs after the internal geometric and thermal properties of the borehole cease to influence the temperature response and before the surface and bottom end effects influence the temperature response. The third transient region occurs when three dimensional effects such as borehole to borehole interaction, surface and bottom end effects influence the temperature response. The borehole transient resistance or gfunction is broken into two zones called the short time step (STS) gfunction (Yavuzturk, 1999) and the long time step (LTS) gfunction. (Eskilson, 1987) The short and the long time step gfunctions relate to the three regions described above in that the short time step gfunction represents region one and two and the long time step gfunction represents region two and three. Thus it is important to note that the short and long time step gfunction can both represent region 2. This allows the two gfunctions to be integrated into one continuous gfunction curve, allowing the borehole transient resistance to be known for small times, such as 0.5 hours, to large times, such as 100 years. Short and long time step gfunctions are discussed in 7 detail in section 1.2.2.3 and 1.2.2.2 respectively. The next two sections, however, will discuss the current literature for steady state and transient borehole modeling. 1.2.1 Steady State Modeling of Boreholes This section discusses borehole resistance since it is an important part of transient analysis. The borehole resistance is the thermal resistance between the fluid and the borehole wall. Figure 13 shows a crosssection of a borehole and a corresponding thermal delta circuit. Figure 13 Crosssection of the Borehole and the Corresponding Thermal ΔCircuit (Hellström, 1991 p.78) f 1 T and f 2 T (°C or °F) represent the fluid temperature in each leg of the Utube and 1 q and 2 q ⎥⎦ ⎤ ⎢⎣ ⎡ hr ⋅ ft or Btu m W the heat flux (heat transfer rate per unit length of borehole) from the circulating fluid. b T represents the average temperature on the borehole wall. As shown in the delta circuit in Figure 13, the thermal resistance between f 1 T and b T is 1 R ⎥⎦ ⎤ ⎢⎣ ⎡ ⋅ ⋅° Btu or hr ft F W mK and the thermal resistance between f 2 T and b T is 2 R 8 ⎥⎦ ⎤ ⎢⎣ ⎡ ⋅ ⋅° Btu or hr ft F W mK . 12 R represents the “short circuit” resistance for heat flow between f 1 T and f 2 T . However, if the fluid temperatures in each leg of the Utube are approximately equal, which occurs at the bottom of the borehole, the resistance 12 R can be neglected in the Δcircuit. The 12 R resistance is often neglected for the entire borehole. This has the effect of decoupling one leg of the Utube from the other, greatly simplifying the system. Figure 14 shows a decoupled borehole system with a circuit diagram defining 1 R and 2 R . GROUT UTUBE GROUND BOREHOLE SYMMETRY LINE RP1 Rf1 Rg1 Rf2 RP2 Rg2 Tf Tb Tb Tf Tf R1 R2 Tb Figure 14 Cross Section of a Borehole with Symmetry Line and the Corresponding Thermal Circuit. In decoupling the borehole, the assumption is made that the grout, pipe, and fluid for each half of the borehole have the same geometry and thermal properties. This 9 assumption means that f 1 f 2 R = R , p1 p2 R = R , and g1 g 2 R = R . Thus 1 R and 2 R , from the circuit in Figure 13, are equal. However, the total resistance of the grout is not typically written in terms of the grout resistance for half of the borehole. The overall grout resistance is instead lumped in one g R term. With these assumptions the borehole resistance circuit shown in Figure 14 can easily be reduced to produce Equation 11. This equation describes the overall borehole resistance. 2 pipe fluid total grout R R R R + = + where, total R = borehole thermal resistance ⎟⎠ ⎞ ⎜⎝ ⎛ W mK or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF grout R = grout thermal resistance ⎟⎠ ⎞ ⎜⎝ ⎛ W mK or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF pipe R = pipe thermal resistance for one tube ⎟⎠ ⎞ ⎜⎝ ⎛ W mK or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF fluid R = fluid thermal resistance for one tube ⎟⎠ ⎞ ⎜⎝ ⎛ W mK or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF (11) The two major contributors to the borehole resistance are the grout and pipe resistance. The fluid resistance contributes typically less than one percent to the overall steady state borehole resistance for turbulent flow. For laminar flow the contribution made by the fluid resistance is much greater and can exceed twenty percent of total R . The pipe resistance can be calculated with Equation 1–2 (Drake and Eckert, 1972). 10 k r r Rpipe 2π ln 1 2 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ = where, pipe R = pipe thermal resistance ⎟⎠ ⎞ ⎜⎝ ⎛ W mK or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF k = Conductivity of the pipe ⎟⎠ ⎞ ⎜⎝ ⎛ mK W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h⋅ ft⋅ F Btu o 2 r = outside diameter (m) or (ft) 1 r = Inside diameter (m) or (ft) (12) The fluid resistance can be calculated using Equation 13 (Drake and Eckert, 1972). r h Rfluid 1 2 1 π = where, fluid R = fluid thermal resistance ⎟⎠ ⎞ ⎜⎝⎛ W mK or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF h = convection coefficient of the fluid ⎟⎠ ⎞ ⎜⎝ ⎛ m K W 2 or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft ⋅ F Btu 2 o 1 r = Utube inside diameter (m) or (ft) (13) The grout resistance can be calculated from the average temperature profile at the borehole wall and the surface of the Utubes with Equation 14 (Hellström, 1991), presuming these temperatures are available. 11 Q R TU tube TBH Wall grout − − − = where, grout R = grout thermal resistance ⎟⎠ ⎞ ⎜⎝ ⎛ W mK or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF Q = Heat flux per unit length of Utube ⎟⎠ ⎞ ⎜⎝ ⎛ m W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft Btu U tube T − = Average temperature at outer surface of Utube (K) or (ºF) BH Wall T − = Average borehole wall temperature (K) or (ºF) (14) The grout resistance is the most complicated component of the borehole resistance. Unlike the pipe and fluid resistance, the apparent grout resistance will change significantly over the first few hours of heat injection or extraction. There are several methods of calculating the grout thermal resistance. The methods that are used to directly calculate the steady state borehole resistance are the Gu and O’Neal (a, 1998) approximate diameter equation, the Paul (1996) method, and the multipole (Bennet and Claesson, 1987) method. Other methods calculate the transient heat transfer between the fluid and surrounding ground, but may be applied to calculate the borehole resistance, include the cylinder source (Ingersoll, 1948) method, and the finite volume method (Yavuzturk, 1999). 1.2.1.1 Line Source Model The line source developed by Kelvin and later solved by Ingersoll and Plass (1948), is the most basic model for calculating heat transfer between a line source and the earth. In this model the borehole geometry is neglected and modeled as a line source or sink of infinite length, surrounded by an infinite homogenous medium. Thus, with 12 respect to modeling a borehole, the line source model neglects the end temperature effects. Equation 15 is the general equation that Ingersoll and Plass (1948) used to model the temperature at any point in an infinite medium from a line source or sink. The medium is assumed to be at a uniform temperature at time zero. ∫ ∞ − Δ = soil x e d k T q β π β β 4 where, t x r soil 4α 2 = (15) (16) ΔT = change in ground temperature at a distance r from the line source (°C) or (°F) q = heat transfer rate per length of line source ⎟⎠ ⎞ ⎜⎝ ⎛ m W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h⋅ ft Btu t = time duration of heat input q (s) r = radius from the line source (m) or (ft) soil α = soil thermal diffusivity ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ s m2 or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ s ft2 soil k = conductivity of the soil ⎟⎠ ⎞ ⎜⎝⎛ mK W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft⋅ F Btu o The integral in Equation 15 can be approximated with Equation 17 for an mth stage of refinement. 13 Σ= ⋅ − = − − − m n n n n n I x x x 1 ! ( ) γ ln( ) ( 1) where, x = Defined by Equation 1–6 γ = 0.5772156649 = Euler’s Constant (17) Equation 17 shows the general form of the line source for the mth stage of refinement. In most references the second stage of refinement is used. This method is only accurate for large times. For a typical borehole this equates to times greater than approximately 10 hours. For small times, less than 10 hours, the GaussLaguerre quadrature approximation, as shown in Equation 18, is given. This approximation uses the fourth order GaussLaguerre quadrature to solve the infinite integral in Equation 15. ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ + + ⋅ + + ⋅ + + ⋅ + = − ⋅ 44 44 43 43 42 42 41 41 ( ) 1 1 1 1 x z w x z w x z w x z I x e x w quad where, 41 w = 0.6031541043 41 z = 0.322547689619 42 w = 0.357418692438 42 z = 1.745761101158 43 w = 0.0388879085150 43 z = 4.536620296921 44 w = 0.00053929470561 44 z = 9.395070912301 (18) The Gauss – Laguerre quadrature approximation shown here should be used for small times, approximately less than 10 hours. Equation 19 is a modification of Equation 15 and shows how to use the line source to model the borehole fluid temperature. Without the borehole resistance ( bh q ⋅ R ) 14 the borehole temperature (T ) would be the temperature at the borehole wall radius and not the fluid temperature. bh ff soil bh soil q R T t I r k T t q + ⋅ + ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ = ⋅ 4π 4α ( ) 2 (19) T = Borehole fluid temperature (°C) or (°F) t = time duration of heat input (s) ff T = far field temperature of the soil (°C) or (°F) q = heat transfer rate per length of line source ⎟⎠ ⎞ ⎜⎝ ⎛ m W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h⋅ ft Btu bh r = radius of the borehole (m) or (ft) bh R = steady state borehole resistance ⎟⎠ ⎞ ⎜⎝ ⎛ W mK or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF soil α = soil thermal diffusivity ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ s m2 or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ s ft2 soil k = conductivity of the soil ⎟⎠ ⎞ ⎜⎝ ⎛ mK W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft⋅ F Btu o Equation 19 differs from Equation 1.5 in that it uses the steady state borehole resistance to model the heat transfer from the borehole wall to the fluid; the line source model is used to model the heat transfer between the borehole wall and the far field. This usage of the line source requires that the steady state borehole resistance is known. Since the steady state resistance is used the line source will have error for short times before the borehole reaches steady state resistance. For most boreholes the error in the steady state borehole resistance is negligible at 2 hours. The line source model is very easy to use and requires relatively few calculations compared to other methods. However, the drawback to this model is that the borehole 15 internal geometry, thermal properties, and the mass of the fluid are not modeled. The resulting inaccuracies will be examined in Chapter 4. 1.2.1.2 GuO’Neal Equivalent Diameter Model The Gu and O’Neal (1998 a) equivalent diameter method is a very simple method of calculating the steady state borehole thermal resistance. It yields a steady state borehole resistance value that is adequate for most simple calculations. This method is represented by an algebraic equation for combining the Utube fluid into one circular region inside the center of the borehole such that the resistance between the equivalent diameter and borehole wall is equal to the steady state borehole resistance of the grout. Equation 1–10 is used to calculate the equivalent diameter. As can be seen the equivalent diameter is based solely on the diameter of the Utube and the center to center distance between the two legs. eq s D = 2D⋅ L s BH D ≤ L ≤ r where, eq D = Equivalent diameter (m) or (ft) BH r = radius of the borehole (m) or (ft) D = diameter of the Utube (m) or (ft) s L = center to center distance between the two legs (m) or (ft) (110) Figure 15 shows three actual configurations and their equivalent diameters. “d” shows the equivalent diameter for configuration “a”; “e” for “b”; and “f” for “c”. 16 (a) (b) (c) (d) (e) (f) Figure 15 Actual Geometry vs Equivalent Diameter Approximation To calculate the grout resistance, Equation 111, which is the general equation for radial heat conduction through a cylinder, should be employed. grout eq bh grout k D D R 2π ln ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = where, grout R = grout thermal resistance ⎟⎠ ⎞ ⎜⎝ ⎛ W mK or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF grout k = conductivity of the grout ⎟⎠ ⎞ ⎜⎝ ⎛ mK W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h⋅ ft⋅ F Btu o bh D = diameter of the borehole (m) or (ft) eq D = equivalent diameter using GuO’Neal’s method (m) or (ft) (111) The steady state resistance of the grout can be used in Equation 1.1 to calculate the overall borehole resistance. 17 1.2.1.3 Paul Model An experimentally and analytically based method for calculating steady state borehole resistance was developed at South Dakota State University by Paul (1996). The Paul method for calculating the steady state borehole resistance was created using both experimental data and a two dimensional finite element program for modeling a borehole cross section. Several different borehole parameters were modeled such as shank spacing, borehole diameter, Utube diameter, grout conductivity, and soil conductivity. The test apparatus used a single layer thick coil of wire wrapped around each side of the Utube to form an electrical resistance heater. This provided a uniform, constant flux heat input for the system. A real borehole will not have uniform flux at the pipe wall. Heat was input until steady state temperature conditions at the borehole wall radius and along the circumference of the Utube were reached. The borehole resistance was then calculated from the temperatures and the flux. A two dimensional finite element model was created using ANSYS, a UNIX based software package, for the purpose of extending the range of borehole diameters and pipe sizes that the steady state borehole resistance could be solved for. The ANSYS cases could be run much faster than the experimental apparatus; this allowed for more cases to be run. Experimental results from the test apparatus and the ANSYS model were compared for validation purposes. From the results, shape factor correlations were created to model the complex geometry of the borehole. Equation 112 is the resulting shape factor equation for calculating the steady state grout resistance. 18 K S R grout grout ⋅ = 1 , 1 0 β β ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ = ⋅ U−tube b d S d where, grout R = Equivalent diameter (m) or (ft) grout K = Conductivity of the grout ⎟⎠ ⎞ ⎜⎝ ⎛ mk W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft⋅ F Btu o S = Shape factor (dimensionless) 0 β and 1 β = Curve fit coefficients (dimensionless) b d = Diameter of the borehole (m) or (ft) U tube d − = Outside diameter of the Utube (m) or (ft) (112) Equation fit coefficients are given by Paul (1996) for four different shank spacings; A0, A1, B and C. The shank spacings are described in Figure 16. Index Spacing Condition A0 S1 = 0 A1 S1=.123 in (.3124 cm) B S1 = S2 C S2 = 0.118 in (.300 cm) S1 S2 Figure 16 Types of shank spacing used in the Paul borehole resistance approximation. The 0 β and 1 β values for the shank spacing described in Figure 16 are given in Table 1 6. Table 11 Paul Curve Fit Parameters used to Calculate the Steady State Grout Resistance A0 A1 B C 0 β 14.450872 20.100377 17.44268 21.90587 1 β 0.8176 0.94467 0.605154 0.3796 R 0.997096 0.992558 0.999673 0.9698754 19 The R value indicates the accuracy of the curve with respect to the experimental or ANSYS model. An R value of 1 indicates a perfect fit. The grout resistance found using this method should be applied within Equation 11 to determine the overall borehole resistance. 1.2.1.4 Cylinder Source Model The cylinder source model, created by Ingersoll and Plass (1948), uses an infinitely long cylinder inside an infinite medium with constant properties and solves the analytical solution of the 2D heat conduction equation. The cylinder source solution for the gfunction and temperature change at the borehole wall can be calculated with Equations 113, 14, and 15. ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ Δ = ⋅ o o soil r G F r k T q , , r 2 t F soil o ⋅ = α ( ) ( ) ( ) ( ( ) ( )) β β β β β β β β π β d J Y r Y J Y r r J r e r G F r F o o o o o ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ + ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ − ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ − = ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ∫ ∞ − ⋅ 2 2 1 2 1 0 1 1 0 0 2 , 1 1 2 where, (113) (114) (115) ΔT = temperature difference between the steady state temperature of the ground and the temperature at the borehole wall (ºC) or (ºF) q = heat flux per unit length of the borehole ⎟⎠ ⎞ ⎜⎝ ⎛ m W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft Btu o F = Fourier number (dimensionless) r = inner cylinder radius (equivalent Utube radius) (m) or (ft) o r = outer cylinder radius (borehole radius) (m) or (ft) 0 J , 1 J , 0 Y , 1 Y = Bessel functions of the zero and first orders t = time (s) 20 The variable “r” is the location at which a temperature is desired from the cylinder source located at o r . G( o F , r / o r ) is a function of time and distance only. To apply the cylinder source equation for modeling the fluid temperature within a borehole, Equation 116 can be used, setting r equal to the equivalent Utube radius and o r equal to the borehole radius. bh ff o o soil q R T r G F r k T t q + ⋅ + ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ( ) = ⋅ , where, (116) T(t) = borehole fluid temperature (ºC) or (ºF) q = heat flux per unit length of the borehole ⎟⎠ ⎞ ⎜⎝ ⎛ m W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft Btu o F = Fourier number (dimensionless) r = inner radius (equivalent Utube radius) (m) or (ft) o r = outer cylinder radius (borehole radius) (m) or (ft) ff T = far field temperature of the soil (ºC) or (ºF) t = time (s) bh R = steady state borehole resistance ⎟⎠ ⎞ ⎜⎝ ⎛ W mK or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF soil k = soil conductivity ⎟⎠ ⎞ ⎜⎝ ⎛ mK W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF The cylinder source can be used to model the steady state borehole resistance using Equation 11. The Utube and fluid resistances can be calculated as shown in section 1.2.1 in Equations 12 and 13. 21 A comparison of borehole resistance calculation methods, including the cylinder source method, is shown in chapter 2. Also a correction factor was created that increases the accuracy of the cylinder source greatly for boreholes with a large shank spacing. As will be shown in chapter 2, even without the correction factor, the cylinder source with superposition technique is a reasonably accurate method for calculating the grout resistance. 1.2.1.5 Multipole The multipole (Bennet, et al. 1987) method is used to model conductive heat flow in and between pipes of differing radius. In the multipole model, the tubes are located inside a homogenous circular region that is inside another homogenous circular region. The multipole method is not constrained to calculating the steady state borehole resistance for a borehole with only one Utube. Furthermore, the tubes do not need to be symmetrical about any axis. This is advantageous since some boreholes have two Utubes. The model is also able to calculate borehole resistance for Utubes that are not equidistant from the center of the borehole. To show the capabilities of the model Figure 17 has been created showing an asymmetric borehole with three pipes. The pipes have temperatures, f 1 T , f 2 T , and f 3 T . 22 Soil Grout Pipe Fluid Tf1 Tf2 Tf3 Ks 1 Kg b 2 b b3 x rb rs 1 r 2 r r 3 bc Tc Figure 17 Example of a 2D System for the Multipole Method. The inner circular region represents the grout and the outer region represents the soil for the borehole system. For calculating borehole resistance, b r can be set to 100 m (328 ft). The inputs to the multipole method are shown in Table 12. 23 Table 12 Variable Input List for the Multipole Method. g K Thermal conductivity in the inner region ⎟⎠ ⎞ ⎜⎝ ⎛ mk W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft⋅ F Btu o s K Thermal conductivity in the outer region ⎟⎠ ⎞ ⎜⎝ ⎛ mk W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft⋅ F Btu o N Number of pipes J Order of multipole b r Radius of the outer region (m) or (ft) s r Radius of the inner region (m) or (ft) c β Thermal resistance coefficient at the outer circle (nondimensional) c T Temperature of the outer region (K) or (ºF) The following are input for each pipe indexed by i i i x , y Location of each innermost pipe (m) or (ft) i r Radius of each pipe (m) or (ft) i β Thermal resistance coefficient for each pipe (nondimensional) fi T Fluid temperature (K) or (ºF) In Table 12 the nondimensional variable i β is used to input the pipe thermal resistances. This is shown in Equation 117. 24 β = 2πRk where, R = Thermal resistance of the pipe ⎟⎠ ⎞ ⎜⎝ ⎛ W Km or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF k = Grout conductivity ⎟⎠ ⎞ ⎜⎝ ⎛ mk W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft⋅ F Btu o β = resistance coefficient (Nondimensional) (117) The general equation that the multipole method solves is the steady state twodimensional heat conduction equation, Equation 118: 0 2 2 2 2 = ∂ ∂ + ∂ ∂ y T x T where, T = Temperature (°C) or (ºF) x = Distance in the x direction (m) or (ft) y = Distance in the y direction (m) or (ft) (118) When solving the differential equation several assumptions are made. The temperature is constant inside of the pipe walls and the fluid convective resistance. The temperature around the outer region is constant. The system is at steady state. Multipoles were created using the line source model. They are called multipoles because for each line source there is a line sink at a mirror point. This can be seen in Equation 119, the temperature equation for a zeroth order multipole, where the line sink n q is at ( , ) n n x y in the first term and the mirror sink of strength n σ ⋅ q is located at the mirror point ( 2 / 2 , 2 / 2 ) n b n n b n x r r y r r in the second term. A zeroth order multipole has one source and one sink. 25 0 ≤ r ≤ rb : (119) ⎪ ⎪ ⎭ ⎪ ⎪ ⎬ ⎫ ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ − + − + ⋅ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − + − = 2 2 2 2 2 2 2 2 ( ) ( ) ln / ( ) ( ) ln 2 ( , ) n b n n b n b n n n b b n r y y r r x x r r r x x y y r k T x y q σ π . k k k k b b + − σ = , −1≤σ ≤1 where, T(x, y) = Temperature (K) or (ºF) n q = Heat flux per unit length ⎟⎠ ⎞ ⎜⎝ ⎛ m W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft Btu b k = Conductivity of the inner region ⎟⎠ ⎞ ⎜⎝ ⎛ mk W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft⋅ F Btu o k = Conductivity of the outer region ⎟⎠⎞ ⎜⎝ ⎛ mk W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft⋅ F Btu o n x = Location of the line source in the x direction (m) or (ft) n y = Location of the line source in the y direction (m) or (ft) b r = Radius of the inner region (m) or (ft) n r = Radius of the line source (m) or (ft) To give a graphical perspective on the location of the source and sink Figure 1.8 is given. As can be seen in Figure 1.8 the sink lies on the same radius line as the source. 26 Sink Source y x Grout Soil Radial line from borehole center Borehole Pipe Figure 18 Source and Sink Location for a Single Pipe Equation 119 is the zeroth order multipole equation; to produce a more accurate solution more sources and sinks can be utilized for each pipe. To do this requires a simplification of Equation 119 by use of polar notation. Writing Equation 119 in polar notation yields Equation 120 where ( ) 0 Re n W is the real component of the zero order multipole. 27 ( ) 0 Re 2 ( , ) n b n W k T x y = q ⋅ π where, T(x, y) = Temperature (K) or (ºF) n q = Heat flux per unit length ⎟⎠ ⎞ ⎜⎝ ⎛ m W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft Btu b k = Conductivity of the inner region ⎟⎠ ⎞ ⎜⎝ ⎛ mk W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft⋅ F Btu o ( ) 0 Re n W = Real component of the zero order multipole (120) For higher order multipoles, derivatives are taken of the n0 W as shown in Equation 121. ( ) ( 1)! 1 j no n j nj W j z W ∂ ∂ ⋅ − = where, nj W = jth order multipole of nth line source j = Order of multipole n z = Location of pipe n in polar coordinates no W = zero order multipole (121) Both the real and imaginary components of nj W satisfy the continuous radial flux boundary condition at b r = r . The constant temperature condition of each of the pipes and the outer radius s r is satisfied using a Fourier series expansion. Using this method the temperature for any point within s r can be found. Equation 122 and 23 show the 28 general solution for the temperature inside the outer cylinder radius for orders of multipoles greater than zero. For a borehole this becomes the borehole wall temperature. ⎥⎦ ⎤ ⎢⎣ ⎡ = + Σ ⋅ +ΣΣ ⋅ ⋅ +Σ ⋅ − ⋅ = = j cj j nj cj c N n j j nj pn N n o n no T T P W P r W P r W 1 1 Re b n n k P q ⋅ = 2π where, (122) (123) T = Temperature at the borehole wall (K) or (ºF) o T = Far field temperature (K) or (ºF) n = counter variable j = Order of multipole N = number of pipe no W = zero order multipole cj W = Multipole of the outer cylinder nj W = jth order Multipole of nth line source pn r = Radius of the innermost pipes (m) or (ft) c r = Radius of the cylinder encircling the innermost pipes (m) or (ft) The final equations, shown in Chapter 8 (Bennet, et al. 1987), are an elaboration of Equation 122. In the paper three equations are presented which must be solved iteratively. In Chapter 11 (Bennet, et al. 1987) Fortran 77 code is conveniently given which solves the equations. 29 The multipole method can produce highly accurate steady state temperature profiles of a borehole and soil. Figure 19 shows the two dimensional steady state temperature inside and around a typical borehole. For this figure the borehole fluid temperature was set to 10 °C (18 °F) above a zero far field temperature and a tenth order multipole was used. This method requires a constant temperature boundary condition inside the fluid convective resistance inside the Utubes. This is a reasonable assumption if the fluid in the Utubes is in the turbulent flow regime. Figure 19 Steady State Temperature Field for a Borehole Heat Exchanger 30 A tenth order multipole produces four or five digits of accuracy. Since the multipole method is very fast on computers above 400 megahertz, a tenth order multipole should be used for most problems. As can be seen in Figure 19, the difference in grout and soil conductivity creates a slope discontinuity at the borehole radius, where the inner circular region representing the grout meets the outer region representing the soil. A typical steady state temperature profile at the borehole wall is shown in Figure 110. Steady State Borehole Wall Temperature 0 1 2 3 4 5 6 7 8 0 50 100 150 200 250 300 350 Distance Along the Circumference (Degrees) Temperature (c) Figure 110 Temperature Change Around the Borehole Circumference To calculate the steady state borehole resistance a temperature is specified for each leg of the Utube. The multipole method is then used to calculate a heat flux out of each Utube and the temperature around the circumference of the borehole radius should be averaged. Since the multipole program solves for temperatures very quickly, using 360 points at the borehole radius is feasible and will produce a high degree of accuracy. 31 Once an average temperature at the borehole radius is found and it can then be used with the flux in Equation 124 to find the steady state borehole resistance. U tube1 U tube2 fluid bh BH q q T T R − − + − = where, BH R = borehole thermal resistance ⎟⎠ ⎞ ⎜⎝ ⎛ W mK or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF fluid T = average temperature at the Utube pipe wall (K) or (ºF) bh T = average temperature at the borehole radius (K) or (ºF) U tube1 q − = heat flux from Utube leg one ⎟⎠ ⎞ ⎜⎝ ⎛ m W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft Btu U tube2 q − = heat flux from Utube leg two ⎟⎠ ⎞ ⎜⎝ ⎛ m W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft Btu (124) 1.2.2 Transient Modeling of Borehole Heat Exchangers Transient heat transfer in boreholes occurs when the heat flux entering the borehole through the fluid does not equal the heat flux leaving the borehole via the borehole wall. Borehole transients have a significant effect on the borehole fluid temperature response after any change in heat extraction or rejection rate. For a step change in the heat flux, the time for which transient effects are significant is determined primarily by the grout thermal conductivity and the borehole geometry such as the shank spacing and radius of the borehole. In general, a small grout thermal conductivity or a large borehole radius will lengthen the transient region. In most actual systems, the heat flux applied to a borehole through the circulating fluid changes continuously throughout operation. Thus borehole transients need to be modeled not only at the beginning of a simulation but throughout the simulation. 32 Several analytical two dimensional models exist and have been used to model boreholes such as the line source (Ingersoll and Plass, 1948) and cylinder source (Ingersoll, 1948). These methods have very limited capability for modeling the internal borehole transients especially for transient heat pulse changes in the first 10 hours. This section, describes the buried electrical cable (Carslaw and Jaeger, 1947) analytical model, and presents the General Elliptical Multiblock Solver (GEMS2D) and Yavuzturk’s (1999) pie sector finite volume method programs. The application of the buried electrical cable model to borehole heat exchangers is covered in Chapter 3. This section also covers Eskilson’s (1987) three dimensional model of boreholes and its coupling with the two dimensional analytical or finite volume methods via borehole resistance. 1.2.2.1. Buried Electrical Cable Model The buried electrical cable (BEC) model (Carslaw and Jaeger 1947) is an analytical model used to describe the heat flow out of a cable buried in the ground. An electrical cable consists of three main parts, a metal core surrounded by insulation and then an outer protective sheath. A diagram of a buried electrical cable is shown in Figure 111 along with a circuit diagram of the system. Implicit in this method are the assumptions that the core and sheath thermal capacities have finite thermal capacities but are perfect conductors and that the insulation has negligible heat capacity, but a fixed thermal resistance. The most significant difference between the buried electrical cable model and other analytical models such as the line source and cylinder source is that this model incorporates the thermal capacity of the sheath and core in calculating the temperature profile of the core. As seen in the circuit diagram in Figure 111 the core and sheath 33 thermal capacities are represented as 1 S and 2 S . The insulation resistance is represented as s R . In the circuit a heat flux can be applied creating a temperature differential between the core and soil. The heat flux is absorbed by the capacitances 1 S and 2 S . INSULATION SHEATH GROUND CORE Rs Tsoil Tcore Tsoil, ff Tsoil, ff S1 S2 Figure 111 Diagram of a Buried Electrical Cable and Circuit The analytical equations, given as Equations 125 and 126, for this system are more complicated and require more computational time than the cylinder source and the line source equations. 34 G(t) k T q soil Δ = (125) du u G t a a e a t u ∫∞ ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ − ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ⋅ Δ − = 3 0 3 2 2 2 1 2 2 ( ) 2 1 α π (126) 1 2 1 2 S r C a b p π ρ = , 2 2 2 2 S r C a b p π ρ = , s soil h = 2π ⋅ R k [[ ( ) ( ) ] [ ( ) ( ) ]2 ] 1 2 0 2 1 2 1 2 2 1 2 0 2 1 2 1 2 Δ = u a + a − hu J (u) − a a − hu J (u) + u a + a − hu Y (u) − a a − hu Y (u) where, t = time (s) b r = outer radius of the sheath (m) or (ft) soil k = Conduc1tivity of the soil ⎟⎠ ⎞ ⎜⎝ ⎛ mK W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft⋅ F Btu o ρ = density of the soil ⎟⎠ ⎞ ⎜⎝ ⎛ m3 kg or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ft3 lbm p C = specific heat of the soil ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ kgK J or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ kg⋅ F Btu o s R = insulation thermal resistance ⎟⎠ ⎞ ⎜⎝ ⎛ W Km or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF 1 S = core thermal capacity ⎟⎠ ⎞ ⎜⎝ ⎛ mK J or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ft⋅ F Btu o 2 S = sheath thermal capacity ⎟⎠ ⎞ ⎜⎝ ⎛ mK J or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ft⋅ F Btu o soil α = soil thermal diffusivity ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ s m2 or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ s ft2 0 J , 1 J , 0 Y , 1 Y = Y and J type Bessel functions of zero and first orders 35 Equation 126 will produce a buried electrical cable gfunction for a particular time. It should be noted that this is the only analytical model presented here which takes into account the thermal mass of the heat generation medium, which in the case of a buried electrical cable is the core. However, there is potential for this model to be modified to model a borehole and account for the fluid mass inside a borehole. The application of the BEC analytical equation in modeling a borehole system is discussed in detail in Chapter 3. 1.2.2.2. Gfunction Model: Long Time Step The long time step (LTS) gfunction (Eskilson 1987) represents the nondimensionalized borehole response for times when threedimensional effects such as borehole to borehole interaction, surface and bottom end effects influence the borehole fluid temperature response. Gfunctions are plotted against the natural log of scaled time where the scaling factor is dependent on the depth of the borehole and the soil thermal diffusivity. As developed by Eskilson (1987), the borehole transient resistance or gfunction can be nondimensionalized with respect to the soil and scaled with respect to the steady state borehole resistance to form a gfunction. Both long and short time step gfunctions can be produced using Equation 127 (Eskilson 1987). In Equation 127 the borehole T term represents the timevarying average temperature at the borehole wall and must be calculated with a numerical or analytical procedure. ground T is the far field temperature and usually remains constant. This gfunction represents the nondimensionalized resistance between the ground and borehole wall. Equation 128 includes the borehole resistance term. 36 ( borehole ground ) b soil s T T Q k H r t g t = − ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ 2π , ( ) borehole ground soil BH b soil s T T k R Q k H r t g t π π 2 2 , + − = ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ where, g = gfunction value (dimensionless) Q = flux per unit length ⎟⎠ ⎞ ⎜⎝ ⎛ m W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft Btu soil k = thermal conductivity of the soil ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ft⋅ F Btu o borehole T = average temperature at the borehole wall (°C) or (°F) ground T = far field temperature of the ground (°C) or (°F) (127) (128) As can be seen from Equation 128, the gfunction has two major parameters s t / t and r H b / . For a specific borehole configuration, the first parameter is the major contributor to the gfunction and the second is a factor that corrects the gfunction according to the borehole radius ( b r ) to depth (H). The r H b / correction factor is relatively minor since it changes the gvalues on the order of one percent or less. The main parameter, in Equation 127, that requires significant calculation time is the average temperature at the borehole wall radius ( borehole T ). Gfunctions are plotted against the natural log of nondimensionalized time. The term s t is called the time scale factor, and can be calculated using Equation 129 (Eskilson 1987). 37 soil s t H 9α 2 = where, s t = time scale factor (s) H = depth of the borehole (m) or (ft) soil α = soil thermal diffusivity ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ s m2 or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ s ft2 (129) The gfunction defined by Equation 127 represents the ground thermal resistance and not the borehole resistance. This is beneficial because it allows a single long time step gfunction to be useful for any borehole geometry and soil conductivity as discussed by Eskilson (1987). The gfunction in Equation 128 is only valid for the specific borehole for which the borehole resistance was calculated. An example of a combined long and short time step gfunction for a single borehole system is shown in Figure 112 and 113. In these figures, the gfunction is plotted against log scale time. Figure 112 was created with Equation 127 and Figure 1 13 with Equation 128. Short and Long Time Step Gfunction for a Single Borehole System Without Borehole Resistance 3 2 1 0 1 2 3 4 5 6 7 16 14 12 10 8 6 4 2 0 2 4 ln(t/ts) gfunction Figure 112 Short and long time step gfunction without borehole resistance Short and Long Time Step Gfunction for a Single Borehole System with Borehole Resistance 0 1 2 3 4 5 6 7 8 9 10 16 14 12 10 8 6 4 2 0 2 4 ln(t/ts) gfunction Figure 113 Short and long time step gfunctions with borehole resistance 38 As can be seen in Figure 113 the gfunction approaches zero at small times. This indicates that as time approaches zero the resistance asymptotically approaches zero due to the steady state temperature profile of the soil. Looking at Figure 112 might give the impression, however, that for small times the resistance is negative but this is an illusion created by subtracting the borehole resistance. When the borehole resistance is added back in, as shown in Figure 113, the resistance approaches zero at short times. At large times, Gfunctions will plateau. This occurs because of borehole end effects. Using Equation 130 an average fluid temperature can be calculated if the gfunction is known. ground b soil s borehole T H r t g t k T Q + ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ = , 2π where, borehole T = average temperature at the borehole wall radius (ºC) or (ºF) ground T = far field temperature of the ground (ºC) or (ºF) g = gfunction value (dimensionless) Q = flux per unit length ⎟⎠ ⎞ ⎜⎝ ⎛ m W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ hr ⋅ ft Btu soil k = thermal conductivity of the soil ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ft⋅ F Btu o (130) There are no published analytical solutions that approximate the gfunctions for multiple borehole systems. This is due to multiple borehole systems dependence on not only the depth of the borehole, but also on the distance between each borehole. The interaction between boreholes is difficult or impossible to analytically model. With the 39 boreholes dissipating different unknown amounts of heat, it is difficult to analytically model the average borehole wall temperature of the system. However, it can be resolved using superposition. In order to create long time step gfunctions for multiple borehole systems Eskilson (1987) created a Fortran 77 program which uses a variable mesh finite difference method with cylindrically symmetric coordinates. The program is described in detail for a single borehole in Eskilson (1987). The program has the ability to input a constant heat flux per unit length of borehole and calculates the resulting temperature at the borehole radius ( borehole T ) for various times. The borehole wall temperature ( borehole T ) is then used in Equation 127 to calculate the LTS gfunction. An example of the type of variable mesh grid that Eskilson used is shown in Figure 114 Borehole (i1, j) (i, j) (i, j1) (i+1, j) (i, j+1) Figure 114 Twodimensional radialaxial mesh for a heat extraction borehole in the ground (Eskilson, 1987) In Figure 114 each cell is a rectangular cross section of a ring. Temperatures are calculated at the center of each cell; however, logarithmic interpolation can be used to 40 find the temperature at other points in the soil. Eskilson gives a detailed analysis of appropriate mesh sizes in the radial and axial direction since the mesh determines the accuracy of the solution. In general, small cells provide good numerical accuracy and the temperatures are valid for smaller times. However, small cells create longer computational times for a computer. Eskilson suggests using no smaller cells than necessary for a particular problem. Although computer technology has improved since 1987, mesh size is still important for large simulations. In the examples Eskilson used, the upper part of the borehole is thermally insulated to a depth of 5 meters (16.4 ft) and the overall borehole depth is 115 m (377 ft). Mesh comparisons for short times of 25 years and long times of 237 and 947 years were conducted. In the end, heuristics were created for determining the appropriate mesh size in the radial and axial direction. To solve the thermal performance for a system with multiple boreholes, Eskilson (1987) used the superposition technique. Two different examples are given (Eskilson, 1987): a 4x4 borehole configuration, with 10 m (33 ft) spacing, and a 12x10 borehole configuration, with 4 m (13 ft) spacing, with the simplest type of loading condition where the heat flux at the borehole wall is constant per unit length of the borehole. To validate the accuracy of the program the line source was used in conjunction with superposition. Eskilson’s program was used to produce the LTS gfunction curves shown in Figure 115. It can intuitively be determined that a tighter borehole field will produce more overall resistance and as the boreholes are spaced farther apart, all multiple borehole systems will approach the single borehole case. This can be seen in Figure 1 15. 41 Long Time Step Gfunction for Different Borehole Spacing 0 10 20 30 40 50 60 70 80 90 100 9 8 7 6 5 4 3 2 1 0 1 2 3 ln(t/ts) gfunction 1.5 m (5ft) 3.0 m (10 ft) 4.6 m (15 ft) 6.1 m (20 ft) 9.1 m (30 ft) 24 m (80 ft) Figure 115 Long Time Step gfunction for a 64 Borehole System in an 8x8 Configuration with Varying Borehole Spacing Since none of the internal properties of the borehole are significant, the long time step gfunction might initially seem simple to solve. Because of three dimensional effects, the gfunction for multiple borehole systems is deceptively complicated to solve. 1.2.2.3. Gfunction Model: Short Time Step The short time step (STS) gfunction describes the transients that occur within the borehole before the borehole reaches steady state conditions. For this transient region, the borehole is modeled as having infinite length since surface and bottom end effects can be neglected. The STS gfunction can be approximated using the line source or the cylinder source as described in section 1.2.1.1 and 1.2.1.4, respectively, to calculate a fluid temperature profile versus time which could be input into Equation 131 to yield a 42 gfunction. Equation 131 calculates a gfunction in terms of the fluid temperature and the steady state borehole resistance. Unlike the long time step gfunction, this gfunction is only valid for the specific borehole internal geometry (shank spacing, borehole radius… etc.) and conductivities that the fluid T and BH R was generated with. ( fluid BH ground ) b soil s T R Q T Q k H r t g t = − − ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ( ) 2 , π where, g = gfunction value (dimensionless) Q = flux per unit length ⎟⎠ ⎞ ⎜⎝ ⎛ m W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft Btu soil k = thermal conductivity of the soil ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ft⋅ F Btu o ground T = far field temperature of the ground (°C) or (°F) fluid T = average temperature of the circulating fluid (°C) or (°F) BH R = borehole resistance ⎟⎠ ⎞ ⎜⎝ ⎛ W mK or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF (131) Another method which is capable of generating greater accuracy than analytical methods is the finite volume model (Patankar, 1980). Two programs that have implemented this model will be discussed. The first is called the General Elliptical Multiblock Solver (GEMS2D) and was developed by Rees (2001). Applications of GEMS2D are reported by Spitler, et al. (1999) and Rees, et al. (2002). This program solves the general convection diffusion equation using a boundary fitted grid. GEMS2D is capable of solving both steady state and transient problems. Boundary fitted grids enable GEMS2D to be applied in solving heat transfer problems with complex geometries such as Utubes within a borehole. Figure 116 shows a GEMS2D boundary fitted grid for half of a borehole, since the geometry is symmetrical. 43 Figure 116 Grid for a cross section of a borehole A complicated grid such as that shown in Figure 116 will require several different blocks to be created and then connected together. Each block is composed of many cells. The cells in each block can then be assigned properties such as conductivity and heat capacity. Different cells within a single block can be assigned different properties. A detailed description of how the GEMS2D program was applied to borehole heat conduction with fluid mass is given in Section 4.1. GEMS2D is capable of calculating the steady state borehole resistance and the transient temperature profile at the borehole wall. The gfunction can be calculated from the average borehole wall temperature ( borehole T ) using Equation 131. GEMS2D is written in the Fortran 90/95 language. A grid generation tool was also written to automate the creation of grids for the GEMS2D simulator. Using a text input file, the grid for convectiondiffusion heat transfer problems can be created with the grid generation tool. In the text file, blocks and boundaries are created and thermal properties of each block are specified. After the grid is created, GEMS2D can then be 44 used to simulate the system. The GEMS2D outputs are given in an output text file. Temperatures at each node at each time increment can be given for transient simulations. The second program was developed by Yavuzturk, et al. (1999) and also uses the two dimensional finite volume model, but with a polar grid. It is specifically developed for modeling the heat flow out of a Utube borehole heat exchanger. Like GEMS2D, Yavuzturk’s program is able to model both the transient and steady state solutions for the temperature field within and around a borehole. To model the geometry of a borehole heat exchanger Yavuzturk uses an approximation for the borehole Utube geometry called the pie sector approximation (Yavuzturk and Spitler, 1999). The piesector approximates the cross section of the Utubes via two “pieshaped” wedges. Figure 117 shows the grid that Yavuzturk’s program creates. Only half the borehole is shown since the system is symmetrical. The pie shaped wedge shown in this figure is representative of one leg of the Utube. It is shown bolded, in the figure, while the actual circular Utube geometry is shown for both legs without bolding. 45 Figure 117 Grid for PieSector Approximation (Yavuzturk, 1999) The grid resolution and pie sector approximation for the Utube geometry is determined by an automated parametric grid generation algorithm and is a function of the borehole and Utube pipe geometry. The algorithm matches the inside perimeter of the circular pipe to the inside perimeter of the pie sector and also creates identical heat flux and resistance conditions near the pipe wall between the circular pipe and the pie sector approximation. The fluid resistance is approximated by adjusting the thermal conductivity of the Utube pipe wall. The total radius of the grid is 3.6 m or (12 ft) so that longer simulation times can be conducted. At this radius, the boundary condition is set to a constant far field temperature. The model was primarily written in the Fortran 77 programming language. Inputs to the model are simple since grid generation is automated. The inputs include shank spacing, Utube diameter, borehole diameter, convection coefficient, the volumetric heat 46 capacity of the soil, grout and Utubes, and the conductivity of the soil, grout and Utubes. These inputs are provided in a text file. The outputs are given in an output text file. Since the Utube geometry is not modeled as accurately as in GEMS2D, the piesector approximation will generally be less accurate than GEMS2D. The benefit of Yavuzturk’s program is that it requires approximately half the simulation time as GEMS2D. The resulting simulation model, discussed in detail in Yavuzturk and Spitler (1999), uses Eskilson’s LTS gfunctions simulation methodology with Yavuzturk’s STS gfunctions simulation methodology. The model has been incorporated into a commercially available GLHE design tool called GLHEPRO (Spitler, 2000). GLHEPRO Version 3 is discussed in detail in section 1.2.3. The simulation model has also been implemented and proved useful in several studies (Yavuzturk and Spitler, 2000, Ramamoorthy, et al. 2001, and Chiasson, 1999). Hybrid GSHP systems use other heat rejection equipment, such as cooling towers, fluid coolers, shallow ponds (Chiasson, et al. 2000; Ramamoorthy, et al. 2001) or pavement heating systems (Chiasson, et al. 2000). For example, several operating and control strategies of a cooling tower hybrid GSHP system are discussed in Yavuzturk and Spitler (2000), and are compared with hourly simulations performed in TRNSYS (SEL 1997). A goal of these studies is to show that hybrid GSHP systems can reduce the size of the GLHE system which in turn can reduce the first cost of the system and the necessary land area. 47 Spitler, et al. (2000) which gives a summary of research and developments in grounds source heat pump systems, design, modeling, and applications for commercial and institutional buildings. Presented in this paper are design methodologies for determining hourly and minutely responses for GLHE designs. 1.2.3 GLHEPRO Version 3 Design Tool GLHEPRO Version 3 (Spitler, 2000) combines a Microsoft windows graphical user interface and a ground loop heat exchanger simulation. The software package developed by Marshall and Spitler is based on the methods developed by Eskilson (1987) at the University of Lund, Sweden. GLHEPRO Version 3 has a library of heat pump performance curves and has the capability of adding user defined heat pump performance curves. GLHEPRO Version 3 also has the flexibility of using SI or English units. GLHEPRO Version 3 uses the long time step gfunctions developed by Eskilson (1987). As discussed earlier, Eskilson (1987) created a 2 dimensional finite difference program that omits the internal borehole properties such as shank spacing, grout and Utubes properties. This was done by assuming a steadystate heat transfer process inside the borehole and modeling the transient process outside the borehole using a finite difference technique, so that the temperature at the borehole wall ( borehole T ) was found. Equation 131 can then be used to create the gfunction. GLHEPRO Version 3 uses data from Eskilson’s finite difference program to model the long time step gfunction for over 250 different borehole configurations. The fluid temperature is found by using the temperature at the borehole wall in conjunction with the borehole steady state resistance. The method used in GLHEPRO Version 3 to calculate the short time step response is the 48 line source (Ingersoll, 1948). The Paul (1996) method was used to calculate the borehole resistance for a specific borehole geometry. GLHEPRO Version 3 received inputs for peak and monthly, heating and cooling loads along with borehole internal geometric and thermal properties and borehole configuration. GLHEPRO Version 3 has the capability of outputting the maximum and minimum monthly fluid temperature entering the heat pump and the energy consumption of the system. GLHEPRO Version 3 also has a sizing mode which requires maximum and minimum limits for the entering fluid temperature to the heat pump. The sizing program will find the minimum depth required for a specific borehole configuration to be within the maximum and minimum user defined fluid temperature limits. 1.3 OBJECTIVES The primary objective is to develop and implement a method whereby engineers can more accurately model peakloaddominant systems without time consuming numerical modeling. Implicit within this main objective are the following specific objectives: 1. Determine an appropriate method for calculating the steady state borehole resistance and implement it in GLHEPRO. 2. Enhance shorttimestep (STS) GLHE simulation methodology to account for thermal mass of the fluid to yield more accurate designs via simulations. 3. Develop an automated method for producing the combined short and long time step gfunction. 4. Evaluate the impact of the more accurate gfunction calculation methodologies on the simulation of GLHE systems. 49 5. Evaluate the impact of the more accurate gfunction methodologies on the simulation of GSHP. 50 2 COMPARISON OF BOREHOLE RESISTANCE CALCULATION METHODS A small change in the steady state borehole resistance has a significant impact on the borehole fluid temperature profile. Since the short time step (STS) gfunction is derived directly from the borehole fluid temperature it includes the borehole resistance. In order to be consistent with the long time step (LTS) gfunction, it is necessary to adjust the STS gfunction to subtract the nondimensional temperature rise due to the borehole resistance. This, in turn, requires accurate knowledge of the borehole resistance. Furthermore, as discussed in Chapter 3, the borehole fluid thermal mass (BFTM) model, which is developed in this thesis for calculating the STS gfunction, requires the steady state borehole resistance to be known. Thus, there was a need for comparing different borehole resistance calculation methods for the purpose of choosing one for the BFTM model. This chapter provides a comparison between different methods for calculating the steady state borehole thermal resistance. Both numerical and analytical methods can be used to determine the steady state resistance of the borehole. The numerical methods require much more computational effort but are generally more accurate than approximate analytical methods such as the Gu and O’Neal (a 1998) equivalent diameter method. The general equation for borehole resistance comes from summing the three resistances (fluid, pipe, and grout) between the fluid inside the Utube and the borehole wall as discussed in section 1.2.1. The fluid resistance is typically calculated with a convection correlation. The pipe resistance is determined as a cylindrical conductive resistance. The grout resistance is more difficult to determine, due to the complex geometry. 51 The methods that are compared in this chapter for calculating the grout resistance are the multipole method, the Paul (1996) method, the Gu and O’Neal (a 1998) approximate diameter method, and the cylinder source method (Ingersoll, 1948, 1954). The two numerical programs that are used to calculate the borehole resistance are GEMS2D (Rees, 2001) and Yavuzturk’s (Yavuzturk and Spitler, 1999) pie sector approximation which both use the finite volume method. 2.1 Borehole Resistance Transient and Steady State For the first few hours of constant heat injection or extraction the borehole resistance is transient. Figure 21 shows how the borehole resistance changes over the first 12 hours after a constant flux is applied. Figure 21 was generated from the average fluid temperature and the average temperature at the borehole wall radius using Equation 21. ( ) Q T T R f borehole total − = where, total R = borehole resistance ⎟⎠ ⎞ ⎜⎝ ⎛ W Km or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF f T = Average fluid temperature (K) or (ºF) borehole T = Average temperature at the borehole wall radius (K) or (ºF) Q = flux per unit length ⎟⎠ ⎞ ⎜⎝ ⎛ m W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft Btu (21) This figure comes from a GEMS2D simulation where the borehole was 11.4 cm (4.5 in) diameter, with a 1.6 cm (0.63 in) shank spacing, standard bentonite grout, a pipe 52 conductivity of 0.39 (W/(m·K)) (0.225 Btu/(h·ft·°F)), and a fluid convection coefficient of 1690 (W /(m2K)) (298 (Btu /(h⋅ ft2 ⋅°F))) As can be seen the borehole resistance is almost constant after about 12 hours. For typical boreholes there is usually less than a 2% difference between the steady state value and the value at 10 hours. This is indicated by Figure 22. Resistance vs Time 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0 1 2 3 4 5 6 7 8 9 10 11 12 Time (hours) Resistance (Km/W) Figure 21 Transient Borehole Resistance Profile vs Time Percent Difference of Transient Resistance with Respect to Steady State Resistance vs Time 0 10 20 30 40 50 60 70 80 90 100 0 1 2 3 4 5 6 7 8 9 10 11 12 Time (hours) Percent Difference (%) Figure 22 Percent Difference in Transient Borehole Resistance with respect to Steady State borehole resistance 53 The rate at which the resistance approaches the steady state value is dependent on the geometry and thermal properties within the borehole. A borehole with a low grout or pipe conductivity requires more time for the borehole to reach steady state. 2.2 Borehole Resistance Calculation from Analytical and Empirical Methods The analytical methods that are compared in this thesis include the multipole method, the Gu and O’Neal approximate diameter method, and the cylinder source method. The Paul method is also included in this section because it is based on curve fits to numerical and analytical data and is not strictly numerically based. The literature review in sections 1.2.1.1 through 1.2.1.4 describes how each method, except for the cylinder source, can be used to model the steady state borehole resistance. An application of the cylinder source for calculating the steady state borehole resistance is described in this section. The Gu and O’Neal method and the Paul method are much simpler to calculate than the multipole and cylinder source methods; however the multipole and cylinder source methods produce more accurate solutions. The Gu and O’Neal method was used exactly as described in section 1.2.1.2 and the Paul method was used exactly as described in section 1.2.1.3. Thus, with regards to these methods, no further explanation is necessary in this chapter. However, because of the complexity of the cylinder source and multipole methods, additional explanation is provided here in addition to what has been previously described in the literature review in sections 1.2.1.4 and 1.2.1.5 respectively. 54 The cylinder source solution can be used to model steady state resistance by using the principal of superposition. To model the grout resistance using Equation 14, the temperature rise from the flux exiting each leg of the Utube can be superimposed to calculate the average Utube outside wall temperature and also the average borehole wall temperature. Implicit in this method is the assumption that the soil conductivity, which is different from grout conductivity, has relatively little influence on the borehole thermal resistance. This method will be explained in greater detail using Figure 23 which shows the cylinder source locations for the calculation of the Utube average outside wall temperature. Figure 23 (a) shows a cylinder source located at (x,y) = (0,0) in an infinite medium of grout. The circle labeled borehole in Figure 23 (a) is shown to indicate the cylinder source location inside the borehole. The cylinder source in Figure 23 (a) is the location where the Utube temperature ( U tube T − ) will be calculated for Equation 14. A circle showing the borehole radius is shown, however the cylinder source method does not make a distinction between conductivities or thermal properties nor does it account for the existence of the other Utube since the medium is infinite. For the purpose of calculating borehole resistance, the conductivity of the soil ( soil k ) in Equation 113 should be replaced with the conductivity of the grout ( grout k ). In Figure 23 (a) the temperature rise at the U–tube radius should be calculated using r = o r in Equation 113. Since the temperature rise will not vary it is unnecessary to calculate several temperatures along the Utube circumference and average them. 55 Figure 23 Cylinder Source Diagram for Calculating the Utube Outside Wall Temperature for use in the Steady State Borehole Resistance Calculation Figure 23 (b) is similar to (a) except it models the cylinder source at the other leg of the Utube. Similar to Figure 23 (a), the cylinder source in (b) is also surrounded by an infinite medium of grout. To show where the temperature rise will be calculated, a circle is drawn showing Utube leg 1, however neither the fluid nor the Utube leg 1 pipe exists in the infinite and continuous medium surrounding the cylinder source. A larger 56 view of the Utubes in Figure 23 (b) is shown in Figure 23 (c) and the geometry of how to calculate the radius (r) for Equation 113 is shown in Figure 23 (d). The temperature rise from the cylinder source at Utube leg 2 should be calculated for several points along the circumference of Utube leg 1. The average temperature should then be calculated. This can be done by calculating “r” with the equation shown in Figure 23 (d) for several different Φ angles. To calculate the overall temperature rise at the Utube radius the average temperature increase from each Utube cylinder source should be superimposed to yield the overall temperature. Equation 25 calculates the resistance between the Utube OD and radius infinity. In a real system this is analogous to calculating the combined grout and soil resistance except the soil properties are the same as grout. Likewise Equation 26 calculates the resistance between borehole OD and radius infinity also using grout properties. Equation 27 finds the resistance of the grout that is located between the Utube and the borehole radius by subtracting the resistance calculated in Equation 26 from that calculated in Equation 25. In a similar manner the average temperature at the borehole radius can be found by averaging the temperature rises created by a cylinder source located at each Utube leg. If the borehole is symmetrical down the middle then the calculation shown in Equation 26 will calculate the resistance between the borehole wall and an infinite radius. 57 r (θ ) d 2 r 2 2r d cos(θ ) BH borehole borehole = + − r (θ ) (d / 2)2 r 2 2r (d / 2)cos(θ ) U tube UT UT = + − − ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ = ⋅ k r G F r k R r r k o o o ( , ,α , ) 1 , , Σ− = − − ⎥⎦ ⎤ ⎢⎣ = ⋅ ⎡ ⋅ + 1 0 ( ( 360), , , ) ( , , , ) 2 1 N U tube u tube UT grout grout UT UT grout grout r k R r r k N R r N R θ θ α α Σ− = ⎥⎦ ⎤ ⎢⎣ = ⋅ ⎡ ⋅ 1 0 1 ( ( 360), , , ) N BH BH UT grout grout r k N R r N R θ θ α grout U tube BH R = R − R − where, (22) (23) (24) (25) (26) (27) grout R = borehole grout resistance ⎟⎠ ⎞ ⎜⎝ ⎛ W mK or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF U tube R − = resistance between Utube OD and radius infinity ⎟⎠⎞ ⎜⎝ ⎛ W mK or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF BH R = resistance between borehole OD and radius infinity ⎟⎠ ⎞ ⎜⎝ ⎛ W mK or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF N = number of points to calculate the resistance around the circumference of the borehole or the Utube (dimensionless) θ = dummy variable used for counting from 0 to N (dimensionless) UT r = radius of the Utube (m) or (in) grout α = thermal diffusivity of the grout ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ s m2 or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ s ft 2 grout k = grout conductivity ⎟⎠ ⎞ ⎜⎝ ⎛ mk W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft⋅ F Btu o R(r, r , , k) o α = general function for resistance as a function of G() from Equation 115 (θ ) U tube r − = radial distance from the center of the cylinder source and the outside diameter of the other leg of the Utube (m) or (in) (θ ) BH r = radial distance from the center of the cylinder source and the outside diameter of the other leg of the Utube (m) or (in) d = the distance between the two centers of the left and right Utube legs (m) or (in) borehole r = borehole radius (m) or (in) 58 As shown in Equation 27, to calculate the grout resistance the resistance between the borehole wall and infinity is subtracted from the resistance between the Utube OD and infinity. Thus, Equations 22 through 27 can be applied to yield a solution for the steady state borehole resistance. In order to use the cylinder source to calculate steady state resistance a Fourier number should be calculated and the number of points to solve for around the borehole and Utube radiuses should be established. The time that was used to calculate the Fourier number was 80,000 hours which causes the solution to converge to five or more digits. The number of points that was chosen for finding the average temperature of both the Utube and borehole wall was 90. As the number of points was increased from 8 to180, the solution converged to five or more digits at 90 points. Other parameters that were chosen were the integration bounds in Equation 115 which are between zero and infinity. It was found that an upper integration bound of 10,000 produces 4 or more significant digits of convergence as compared to an integration bound of 100,000 which gives more than 8. The multipole resistance was found using a modified version of the Fortran 77 source code given in Bennet and Claesson (1987). Within the multipole method, the borehole resistance is found by establishing a temperature at the Utube wall and then calculating a heat flux and a temperature profile around the circumference of the borehole wall. The temperature at the borehole was calculated by taking an average of 180 points along the circumference of the borehole wall. Averaging 180 points versus averaging 360 points produced a temperature difference of less than 0.00001 °C (0.000018 °F) difference. The resistance can be calculated by using Equation 11. 59 2.3 Borehole Resistance Calculation using Numerical Methods As described in section 1.2.2.3, GEMS2D closely approximates the borehole geometry using a boundary fitted grid, whereas Yavuzturk approximates the borehole geometry using a pie shaped wedge in a parametric grid to represent the Utubes. In this chapter, GEMS2D is used as a standard for the borehole resistance calculation since it correlates very closely with another highly accurate method, the multipole method. This will be shown in section 2.5. GEMS2D has also proven to be a very accurate two dimensional finite volume program for other simulations. The disadvantage of GEMS2D is that it is approximately half as fast as Yavuzturk’s finite volume model program. In both GEMS2D and Yavuzturk’s pie sector approximation the average borehole wall temperature was subtracted from the given fluid temperature. This is shown in Equation 21. With constant flux and large times, typically greater than 10 hours, Equation 21 produces the steady state borehole resistance. A comparison of the steady state borehole resistances that GEMS2D and Yavuzturk’s pie sector approximation produce is shown in section 2.4. 2.4 Numerical Methods: Comparison between GEMS2D and the PieSector Approximation for Calculating Steady State Resistance Table 21 shows the baseline borehole system configuration that was used for the comparison. By varying individual parameters, this borehole configuration produced Table 22 which shows the steady state borehole resistance for both GEMS2D and the pie sector approximation. The conductivities that were varied in this comparison are the soil, grout and pipe conductivity. 60 The conductivities that were shown cover most typical borehole configurations. Also the borehole diameters that were chosen are 11.4 cm (4.5 in), 15.2 cm (6 in) and 19.1 cm (7.5 in) which also cover typical borehole configurations. As the borehole diameter was changed the shank spacing was held constant at 0.16 cm (0.067 in). The final parameter that was varied was the fluid flow rate at 0.000189 m3/s (3 gpm), 0.000379 m3/s (6 gpm), and 0.000568 m3/s (9 gpm). This also covers the range of most boreholes. Table 21 Borehole Properties (Base Case) Borehole System Table English Units SI Units Diameter 4.5 (in) 114.3 (mm) Shank Spacing 0.067 (in) 1.7 (mm) Utube OD 1.05 (in) 26.67 (mm) Utube ID 0.824 (in) 20.93 (mm) soil k 1.2 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft⋅ F Btu o 2.077 ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W o grout k 0.4 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft⋅ F Btu o 0.692 ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W o pipe k 0.8 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft⋅ F Btu o 1.38 ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W o fluid k 0.8 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft⋅ F Btu o 1.38 ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W o Fluid Volumetric Heat 62.4 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ft ⋅ F Btu 3 o 4.18 ⎟⎠ ⎞ ⎜⎝ ⎛ m ⋅K MJ 3 Flow Rate 3 (gpm) 0.000189 (m3 / s) 61 Table 22 Borehole Resistance Comparison between GEMS2D and the pie sector approximation Varied Input Borehole Resistance % Input s English SI Pie Sector GEMS2D Diff. (oF / h⋅ Btu) (K /W) (oF / h⋅Btu) (K /W) % soil k 0.8 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h⋅ ft ⋅°F Btu 1.38 ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W o 0.3615 0.685 0.3603 0.683 0.333 soil k 1.2 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft ⋅°F Btu 2.07 ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W o 0.3608 0.684 0.3588 0.680 0.559 soil k 1.6 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft ⋅°F Btu 2.77 ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W o 0.3604 0.683 0.3580 0.679 0.671 grout k 0.2 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft ⋅°F Btu 0.346 ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W o 0.6778 1.28 0.6481 1.23 4.48 grout k 0.4 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h⋅ ft ⋅°F Btu 0.692 ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W o 0.3608 0.684 0.3588 0.680 0.559 grout k 0.8 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h⋅ ft ⋅°F Btu 1.38 ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W o 0.1982 0.376 0.2143 0.406 7.79 pipe k 0.4 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h⋅ ft ⋅°F Btu .692 ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W o 0.3898 0.739 0.4284 0.812 9.43 pipe k 0.8 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h⋅ ft ⋅°F Btu 1.38 ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W o 0.3608 0.684 0.3588 0.680 0.559 pipe k 1.2 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h⋅ ft ⋅°F Btu 2.07 ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W o 0.3494 0.662 0.3329 0.631 4.84 Dia. 4.5 (in) 11.4 (cm) 0.3608 0.684 0.3588 0.680 0.559 Dia. 6 (in) 15.2 (cm) 0.4259 0.807 0.3138 0.595 30.3 Dia. 7.5 (in) 19.1 (cm) 0.4758 0.902 0.2785 0.528 52.3 Flow Rate 3 (gpm) 0.000189 (m3 / s) 0.3608 0.684 0.3588 0.680 0.559 Flow Rate 6 (gpm) 0.000379 (m3 / s) 0.3585 0.680 0.3570 0.677 0.413 Flow Rate 9 (gpm) 0.000568 (m3 / s) 0.3576 0.678 0.3563 0.675 0.358 Several observations can be made between the steady state resistance obtained from GEMS2D and the pie sector approximation in Table 22. The pie sector approximation deviates from the GEMS2D solution when the grout geometric properties are changed. By changing the diameter of the borehole without changing the shank spacing or the Utube diameter the grout geometry is being changed. This test measures 62 how accurately the automated grid generation algorithm approximates the actual geometry of the borehole with the pie sector approximation. The resistance is overestimated by more than 50% in the 19.1 cm (7.5 in) diameter case. When the conductivity of the pipe or grout are changed from the standard of 1.38 and 0.692 ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W o (0.8 and 0.4 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h⋅ ft ⋅°F Btu ) respectively the relative percent difference increases substantially. This error is accounted for because, when the conductive properties of the grout and pipe are increased or decreased, the effects of the geometric differences between the two programs are magnified. Changing the conductivity of the soil does not appreciably change the percent difference between the two programs. This is because soil conductivity has a second order effect on borehole resistance since it is outside the borehole and both programs accurately represent the circular geometry at the borehole wall radius. Both GEMS2D and the pie sector approximation would be poor choices as the resistance calculator for use with the BFTM model since they are very slow. The pie sector approximation requires on the order of thirty minutes and GEMS2D requires an hour on a 450 MHZ computer. Therefore analytical methods need to be compared to arrive at a reasonable solution. The pie sector approximation will not be considered further since GEMS2D is the more accurate finite volume model program for predicting borehole resistance as shown in section 2.5. 2.5 Comparison of Methods for Calculating Steady State Borehole Resistance The steady state borehole resistance calculation methods that are compared in this chapter include the Paul (1996) method, the Gu and O’Neal (a 1998) approximate 63 diameter method, cylinder source, and the multipole methods. The GEMS2D solution is given as a general comparison. The data for the baseline borehole used in this study are given in Table 23. Two different grout types were chosen, standard bentonite and thermally enhanced grout. This will give an understanding for how grout conductivity affects the different borehole resistance calculation methods. For each grout type, resistances were calculated for three different borehole diameters 7.6 cm (3 in), 11.4 cm (4.5 in), and 15.2 cm (6 in). The 7.6 cm (3 in) diameter case is an unrealistic borehole configuration; however it is useful for testing the capabilities of the models. For each borehole diameter, resistances were calculated for four different Utube shank spacings ranging from 3.2 mm (0.125 in) from the outside wall of each Utube, to where both Utubes are touching the borehole outside wall. The parameters that were not varied include the Utube diameter, Utube thermal properties, soil thermal properties, and the circulating fluid’s convection coefficient. Table 23 Base Line Borehole Properties Borehole Diameter Pipe  1" SDR11 D = 114.30 mm 4.5 in I.D. = 27.4 mm 1.08 in Grout  Standard Bentonite O.D. = 33.4 mm 1.31 in K = 0.75 ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W o .433 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft ⋅°F Btu K = 0.390 ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W o 0.225 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft ⋅°F Btu p C ρ = 3.90 ⎟⎠ ⎞ ⎜⎝ ⎛ m K MJ 3 58.2 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ft ⋅°F Btu 3 p C ρ = 1.77 ⎟⎠ ⎞ ⎜⎝ ⎛ m K MJ 3 26.4 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ft ⋅°F Btu 3 Soil  Typical Properties Spacing K = 2.50 ⎟⎠⎞ ⎜⎝ ⎛ m⋅ K W o 1.44 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h⋅ ft ⋅°F Btu A S1= 3.18 mm 1/8 in p C ρ = 2.50 ⎟⎠ ⎞ ⎜⎝⎛ m K MJ 3 37.3 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ft ⋅°F Btu 3 B S1= S2 S2 Fluid Convection Coefficient C3 S2= 3.00 mm 0.12 in H = 1690 ⎟⎠ ⎞ ⎜⎝ ⎛ m K W 2 298 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft ⋅°F Btu 2 C S2= 0 mm 0 in S1 S2 64 Table 24 gives the resistances that were calculated using the given borehole properties with the different methods. Table 25 shows the percent error of the steady state borehole resistance with respect to the GEMS2D calculated borehole resistance. Since GEMS2D is not capable of calculating the borehole resistance for the C spacing case, the multipole solution was used for the error calculation. The general methods for calculating the borehole resistance were shown in the literature review for the cylinder source, the multipole, Gu and O’Neal, and the Paul methods. The cylinder source column shown in Table 24 shows data using the cylinder source method described in section 2.2. The multipole data shown in Table 24 uses the tenth order multipole solution. 65 Table 24 Steady State Borehole Resistance Comparison N/A signifies that the method was not suitable for calculating the borehole resistance for this case Borehole Diameter Utube Spacing grout k Paul Gu and O'Neal Cylinder Source Multipole GEMS2D (mm) (in) (W / Km) ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅° Btu hr ft F (Km/W) ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅° Btu hr ft F (Km/W) ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅° Btu hr ft F (Km/W) ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅° Btu hr ft F (Km/W) ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅° Btu hr ft F (Km/W) ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅° Btu hr ft F 76.2 (3) A 0.75 (0.4333) 0.188 (0.326) 0.136 (0.235) 0.1337 (0.2314) 0.1213 (0.2099) 0.1211 (0.2095) 76.2 (3) B 0.75 (0.4333) 0.170 (0.294) 0.136 (0.235) 0.1338 (0.2316) 0.1214 (0.2101) 0.1213 (0.2099) 76.2 (3) C3 0.75 (0.4333) N/A 0.135 (0.233) 0.1331 (0.2303) 0.1206 (0.2087) 0.1204 (0.2084) 76.2 (3) C 0.75 (0.4333) 0.127 (0.220) 0.119 (0.206) 0.1170 (0.2025) 0.1025 (0.1774) N/A 114.3 (4.5) A 0.75 (0.4333) 0.256 (0.443) 0.222 (0.384) 0.2197 (0.3803) 0.2119 (0.3668) 0.2116 (0.3663) 114.3 (4.5) B 0.75 (0.4333) 0.205 (0.354) 0.190 (0.329) 0.1882 (0.3258) 0.1823 (0.3155) 0.1822 (0.3154) 114.3 (4.5) C3 0.75 (0.4333) N/A 0.146 (0.252) 0.1437 (0.2487) 0.1288 (0.2230) 0.1288 (0.2230) 114.3 (4.5) C 0.75 (0.4333) 0.141 (0.244) 0.137 (0.238) 0.1355 (0.2346) 0.1149 (0.1989) N/A 152.4 (6) A 0.75 (0.4333) 0.322 (0.557) 0.283 (0.489) 0.2807 (0.4859) 0.2737 (0.4737) 0.2734 (0.4733) 152.4 (6) B 0.75 (0.4333) 0.235 (0.407) 0.227 (0.392) 0.2248 (0.3892) 0.2216 (0.3836) 0.2216 (0.3836) 152.4 (6) C3 0.75 (0.4333) N/A 0.163 (0.282) 0.1611 (0.2788) 0.1386 (0.2399) 0.1387 (0.2401) 152.4 (6) C 0.75 (0.4333) 0.152 (0.263) 0.157 (0.273) 0.1556 (0.2694) 0.1260 (0.2182) N/A 76.2 (3) A 1.5 (0.8666) 0.116 (0.201) 0.0896 (0.155) 0.08779 (0.1520) 0.08796 (0.1523) 0.08774 (0.1519) 76.2 (3) B 1.5 (0.8666) 0.107 (0.185) 0.0896 (0.155) 0.08786 (0.1521) 0.08803 (0.1524) 0.08788 (0.1521) 76.2 (3) C3 1.5 (0.8666) N/A 0.0893 (0.155) 0.08743 (0.1513) 0.08763 (0.1517) 0.08740 (0.1513) 76.2 (3) C 1.5 (0.8666) 0.0853 (0.148) 0.0813 (0.141) 0.07947 (0.1376) 0.07899 (0.1367) N/A 114.3 (4.5) A 1.5 (0.8666) 0.150 (0.259) 0.133 (0.230) 0.1308 (0.2264) 0.1317 (0.2280) 0.1315 (0.2276) 114.3 (4.5) B 1.5 (0.8666) 0.124 (0.215) 0.117 (0.202) 0.1151 (0.1992) 0.1158 (0.2005) 0.1157 (0.2003) 114.3 (4.5) C3 1.5 (0.8666) N/A 0.0946 (0.164) 0.09279 (0.1606) 0.09149 (0.1584) 0.09144 (0.1583) 114.3 (4.5) C 1.5 (0.8666) 0.0922 (0.160) 0.0905 (0.157) 0.08871 (0.1536) 0.08627 (0.1493) N/A 152.4 (6) A 1.5 (0.8666) 0.183 (0.316) 0.163 (0.282) 0.1613 (0.2793) 0.1624 (0.2811) 0.1621 (0.2806) 152.4 (6) B 1.5 (0.8666) 0.139 (0.241) 0.135 (0.234) 0.1334 (0.2309) 0.1345 (0.2328) 0.1344 (0.2327) 152.4 (6) C3 1.5 (0.8666) N/A 0.103 (0.179) 0.1015 (0.1757) 0.09828 (0.1701) 0.09833 (0.1702) 152.4 (6) C 1.5 (0.8666) 0.0978 (0.169) 0.101 (0.174) 0.09876 (0.1710) 0.09413 (0.1629) N/A 66 Tables 24 and 5 show that the tenth order multipole and GEMS2D correlate very closely yielding a maximum difference of 0.26 percent. In addition, the multipole method is very fast with a computer compared to GEMS2D. It takes less than a second to calculate on a 450 MHz computer whereas the GEMS2D program might require half an hour, depending on the grid, on the same computer. Table 25 Percent Error of Borehole Resistance N/A signifies that the method was not suitable for calculating the borehole resistance for this case Bore Diameter Spacing Kgrout Paul Gu and O'Neal Cylinder Source Multipole mm (in) (W / Km) ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅° Btu hr ft F 76.2 (3) A 0.75 (0.4333) 55.5 11.9 10.4 0.17 76.2 (3) B 0.75 (0.4333) 39.9 11.8 10.4 0.11 76.2 (3) C3 0.75 (0.4333) N/A 12.0 10.5 0.17 76.2 (3) C 0.75 (0.4333) 23.8 15.9 14.1 N/A 114.3 (4.5) A 0.75 (0.4333) 20.8 4.69 3.83 0.12 114.3 (4.5) B 0.75 (0.4333) 12.3 4.28 3.29 0.03 114.3 (4.5) C3 0.75 (0.4333) N/A 12.9 11.5 0.01 114.3 (4.5) C 0.75 (0.4333) 22.5 19.5 17.9 N/A 152.4 (6) A 0.75 (0.4333) 17.8 3.35 2.67 0.09 152.4 (6) B 0.75 (0.4333) 6.14 2.31 1.46 0.01 152.4 (6) C3 0.75 (0.4333) N/A 17.5 16.1 0.10 152.4 (6) C 0.75 (0.4333) 20.6 24.9 23.5 N/A 76.2 (3) A 1.5 (0.8666) 32.2 2.10 0.05 0.25 76.2 (3) B 1.5 (0.8666) 21.3 2.00 0.02 0.17 76.2 (3) C3 1.5 (0.8666) N/A 2.13 0.04 0.26 76.2 (3) C 1.5 (0.8666) 7.96 2.86 0.60 N/A 114.3 (4.5) A 1.5 (0.8666) 13.9 0.86 0.50 0.20 114.3 (4.5) B 1.5 (0.8666) 7.25 0.95 0.59 0.07 114.3 (4.5) C3 1.5 (0.8666) N/A 3.43 1.48 0.05 114.3 (4.5) C 1.5 (0.8666) 6.89 4.89 2.82 N/A 152.4 (6) A 1.5 (0.8666) 12.8 0.64 0.48 0.17 152.4 (6) B 1.5 (0.8666) 3.72 0.57 0.77 0.02 152.4 (6) C3 1.5 (0.8666) N/A 5.03 3.22 0.05 152.4 (6) C 1.5 (0.8666) 3.91 6.81 4.92 N/A Compared to the tenth order multipole method, the Gu and O’Neal approximate diameter method, the cylinder source methods as well as the Paul method typically have greatly reduced accuracy. For all of the models in Table 25, the largest errors occur for 67 the 7.6 cm (3 in) diameter Utube case for both thermally enhanced and non thermally enhanced grout. The Paul method, the Gu and O’Neal method, and the cylinder source method all tend to over predict borehole resistance. The cylinder source is in some cases higher and some cases lower than the actual resistance as can be seen in Table 24. For both the Gu and O’Neal method and the cylinder source method, as the shank spacing increases from “A” (narrowly spaced Utube) to “C” (widely spaced Utube) the error increases substantially. For the Gu and O’Neal method, with standard grout, 11.4 cm (4.5 in) diameter with the “A” and “C” shank spacing, borehole errors were 4.7% and 19.5% respectively. For the same condition the cylinder source solution produced errors of 3.8% for the “A” shank spacing and 17.9% for the “C” shank spacing. This increase in error stems from the Gu and O’Neal method and the cylinder source method not taking into account the soil conductivity. As the Utubes move from very close together to very far apart the impact of soil conductivity on the borehole resistance increases. Thus, as would be expected for the thermally enhanced grout cases, the errors have all substantially decreased for both the Gu and O’Neal and the cylinder source methods, due to the grout and the soil conductivities being closer together. As stated earlier, the data in Table 25 shows that the Gu and ONeal method has an increase in error as the shank spacing increases. Thus, since the resistance is a direct result of the equivalent diameter (Equation 111), the data shows that the equivalent diameter calculation is less accurate for large shank spacings versus small shank spacings. The Paul method performed poorly in comparison to the other methods. In most cases the error produced by the Paul method was several times that of the other methods 68 shown in Table 25. Also, the error fluctuates differently with shank spacing than the Gu and O’Neal method and the cylinder source model. As mentioned in section 1.2.1.3 the experimental model had uniform heat flux around the Utubes which is not the case in a real system. This is the cause of some of the error in the Paul method however it probably does not account for all of the error in the 76.2 mm (3 in) diameter cases shown in Table 2.5. 2.6 Borehole Resistance and Merging of the Short and Long Time Step GFunction The steady state borehole resistance parameter is used to separate the long time step gfunction from specific borehole geometries making a single long time step gfunction valid for any specific borehole geometry. This is accomplished by the total R term in Equation 28. ( ) Q k T R Q T g soil f total ground − ⋅ − = 2π ( ) where, g = gfunction (nondimensionalized) total R = borehole resistance ⎟⎠ ⎞ ⎜⎝ ⎛ W Km or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ⋅ ⋅ Btu h ft oF f T = average fluid temperature (K) or (ºF) ground T = steady state ground temperature (K) or (ºF) Q = flux per unit length ⎟⎠ ⎞ ⎜⎝ ⎛ m W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ h ⋅ ft Btu soil k = soil conductivity ⎟⎠ ⎞ ⎜⎝ ⎛ m⋅ K W or ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ft⋅ F Btu o (28) 69 The gfunction in Equation 28 is only a representation of the thermal resistance of the ground. Before the long time step gfunction can be used to calculate the fluid temperature the borehole resistance must be calculated using specific borehole parameters and then added to the thermal resistance of the ground. If the resistance calculation is not accurate then the long and short time step gfunctions will merge poorly. Figure 21 shows a short time step gfunction calculated with the line source for the borehole with properties shown in Table 21 with a 11.4 cm (4.5 in) diameter and B shank spacing. The steady state borehole resistance is shown in Table 24. In Figure 21 the long time step gfunction is for a single borehole. As can be seen in Figure 21 three different curves have been created for the longtime step gfunction using three different methods for calculating borehole resistance. The “LTS: Generalized” curve is the long time step gfunction without the borehole resistance. 70 Line Source Short Time Step Gfunction Compared to Long Time Step Gfunction Translated Using Different Borehole Resistance Calculation Methods 1 0 1 2 3 4 5 6 7 8 14 13 12 11 10 9 8 7 6 5 4 3 2 1 0 1 2 3 Log Time (ln(t/ts)) gfunction STS: Line Source LTS: GEMS2D and Multipole LTS: GuO'Neal LTS: Paul LTS: Generalized Figure 24 Line Source STS Gfunction Compared to LTS Gfunction Using Different Borehole Resistance Calculation Methods for a Single Borehole System As can be seen in Figure 23 the LTS and STS gfunction merges well using the resistance calculated with either the GEMS2D or Multipole resistance methods. Also the LTS gfunction using the Gu and O’Neal or the Paul methods matches less well with the STS gfunction. As shown in Table 25, the errors in the borehole resistances are 12.3% for the Paul method and 4.3% for the Gu and O’Neal method. The percent errors shown for this particular case in Table 25 are not the greatest errors. For some cases the merging between the long and short time step gfunctions will be even worse using the GuO’Neal and the Paul methods. 71 2.7 Conclusion As discussed in the literature review the long and short time step gfunctions are produced using different methods. The long time step gfunction is produced using superposition with data from a two dimensional radialaxial finite difference model. The short time step gfunction is produced with a two dimensional analytical or experimental model of the cross section of the borehole. Before the gfunction can be used in a simulation, consistency must be checked between the two methods that produce the short and long time step gfunctions and borehole resistance. If the short and long time step gfunction do not merge well together this is evidence of a problem with the borehole resistance calculation or with the short or long time step gfunction itself. This study shows that since the Paul method, for most geometries, does not accurately calculate the borehole resistance and therefore does not ensure a good merge of the long and short time step gfunction, it should not be used in simulations. The Gu O’Neal method is superior to the Paul method and might be suitable in a simulation when a very simple method is needed. The user should be aware of the errors involved with this simple calculation as shown in Table 25. Of the methods that are compared in this chapter, the multipole method is the best analytical method for the purpose of merging the long and short time step gfunction. Also, since the borehole resistance for most simulations will only be computed once, for a given simulation, it is not necessary for the resistance calculator to be exceptionally fast. However using the finite volume methods such as the pie sector approximation or GEMS2D which require fifteen minutes and 30 minutes, respectively, on a 1.4 Ghz computer is not practical. Since the multipole method requires less then a second to calculate on a 450 Mhz Pentium II and attains a 72 very good correlation with the GEMS2D model it is a very good choice for the borehole resistance calculator. 73 3 SHORT TIME STEP GFUNCTION CREATION AND THE BOREHOLE FLUID THERMAL MASS MODEL (BFTM) The short time step gfunction can be generated by any program or equation that is capable of approximating a transient borehole fluid temperature profile over time. The simplest and fastest method for use in a comp 



A 

B 

C 

D 

E 

F 

I 

J 

K 

L 

O 

P 

R 

S 

T 

U 

V 

W 


