

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


DEVELOPMENT, VERIFICATION, AND IMPLEMENTATION OF A HORIZONTAL BURIED PIPE GROUND HEAT TRANSFER MODEL IN ENERGYPLUS By EDWIN LEE Bachelor of Science Oklahoma State University Stillwater, Oklahoma 2006 Submitted to the Faculty of the Graduate College of Oklahoma State University in partial fulfillment of the requirements for the Degree of MASTER OF SCIENCE May, 2008 COPYRIGHT c By EDWIN LEE May, 2008 DEVELOPMENT, VERIFICATION, AND IMPLEMENTATION OF A HORIZONTAL BURIED PIPE GROUND HEAT TRANSFER MODEL IN ENERGYPLUS Thesis Approved: Dr. Daniel Fisher Thesis Advisor Dr. Jeffrey Spitler Dr. David Lilley Dr. A. Gordon Emslie Dean of the Graduate College iii ACKNOWLEDGMENTS There are many people who have helped me along the way. Some have directly helped me with my studies or by providing a bit of humor to ease stressful times. Some have provided inspiration to keep moving in the right direction. It is not conceivable that I could write such a brief summary that would express all the gratitude that I hold for all the ways I have been helped, but I’ll give it the old college try. My utmost praise and worship go to our Lord God Almighty. Without him, there would be nothing, and through him my life has been realized. My mother and father raised me to always be grateful for what I have, and always remember that I have brothers in this world, and guardian angels in the next world watching over me. They have always supported my ventures, even when they seemed sudden and unclear. My brother and sister have always shown me dedication which is unmovable. Throughout tough times, they persevered and created a wonderful family that I am very proud of. They show me understanding when my life seems so busy, and they remind me that there is always time for family. My nieces are seemingly proof that energy is not conserved. They are full of jokes, laughter, smiles, and excitement twentyfour hours a day. Every time I see them, I remember that there is so much in this world to be thankful for, and to never take time for granted. My grandparents have been great rocks which I have clung to. They have always shown courage, and they have always made an effort to make sure I made the right decisions in life. Of my four grandparents, two have passed on, and I miss them very iv much, but I always keep their memories in my heart, and I know they watch over me. My advisors throughout my collegiate career have had drastic influences in my life. Dr. Cremaschi gave me my first work opportunity during my master’s degree, and I learned a great deal from him and my experiences. Dr. Spitler has always been a tough professor, but his lessons are all worth learning, and I have gained invaluable tools throughout my work with him. My graduate advisor, Dr. Fisher, was also my undergraduate advisor. He was the person who opened my eyes to the thermal sciences and created an excitement within me to pursue my master’s degree in this field. His enthusiasm for this field spreads everywhere he goes, and I thank him for motivating me at many times to keep going forward. My master’s degree was not all about hard work, I enjoyed making some great friends throughout this time as well. Sankar, Xiaowei, Steve, Jason, Sam, Michael, Chris and Yang have been there to either play frisbee, soccer, printer tossing, or just to enjoy some dinner with friends. I have always looked forward to their camaraderie. The most influential person in my life is also the love of my life. My wife has been with me for over four years at the time of this writing, and motivated me to continue working no matter how tough it seemed. I have needed her kind words and unwavering love throughout these past years. I have experienced a deep recommitment to Christ since meeting her, and I couldn’t be more proud to call her my wife. v TABLE OF CONTENTS Chapter Page 1 Introduction 1 2 Literature Review 4 2.1 SemiAnalytic Solutions . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.2 Numerical Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.3 Response FactorModels . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.4 ExperimentalWork . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.5 Earth Air Heat Exchanger Analysis . . . . . . . . . . . . . . . . . . . 14 2.6 District Heating . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.7 Soil Property Prediction . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.8 Transport Delay . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.9 Literature Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3 Fluent Reference Model 23 3.1 FluentModel Description . . . . . . . . . . . . . . . . . . . . . . . . 23 3.2 FluentModel Grid Refinement . . . . . . . . . . . . . . . . . . . . . . 24 4 EarthTubeModel 28 4.1 Earth TubeModel Description . . . . . . . . . . . . . . . . . . . . . . 28 4.2 Model Comparisons: VBA Earth Tube vs. Fluent . . . . . . . . . . . 33 4.3 Model Comparison Results: VBA Earth Tube and New Fluent Earth Tube . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 4.4 Earth Tube Conversion to PipeModel . . . . . . . . . . . . . . . . . 40 vi 4.5 Earth TubeModel Guidance . . . . . . . . . . . . . . . . . . . . . . . 45 4.6 Earth TubeModel Conclusions . . . . . . . . . . . . . . . . . . . . . 48 5 MeiModel 49 5.1 Model Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 5.2 Modifications toMeiModel . . . . . . . . . . . . . . . . . . . . . . . 62 5.3 Model Testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 5.3.1 Fully ConstantModel Testing . . . . . . . . . . . . . . . . . . 63 5.3.2 ConstantWall Temperature Testing . . . . . . . . . . . . . . . 63 5.3.3 Initial Condition Independence Testing . . . . . . . . . . . . . 66 5.3.4 High Flow Rate Testing (2DModel Simulation) . . . . . . . . 68 5.4 MeiModel Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . 70 6 PiechowskiModel 73 6.1 Model Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 6.2 Simplified Piechowski Model . . . . . . . . . . . . . . . . . . . . . . . 83 6.3 Analysis of Piechowski Model Grid Dependence Parameters . . . . . . 85 6.4 InitialModel Testing . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 6.4.1 Code Integrity Test . . . . . . . . . . . . . . . . . . . . . . . . 90 6.4.2 Steady State Test . . . . . . . . . . . . . . . . . . . . . . . . . 90 6.4.3 Initial Condition Independence Testing . . . . . . . . . . . . . 94 6.4.4 High Flow Rate Test . . . . . . . . . . . . . . . . . . . . . . . 94 6.5 Comparison to VerificationModel . . . . . . . . . . . . . . . . . . . . 95 6.6 Piechowski Model Conclusions . . . . . . . . . . . . . . . . . . . . . . 102 7 Model Selection 103 7.1 IndividualModel Results Summary . . . . . . . . . . . . . . . . . . . 103 7.1.1 Transfer FunctionMethods . . . . . . . . . . . . . . . . . . . . 103 7.1.2 Earth Tube PipeModel . . . . . . . . . . . . . . . . . . . . . 104 vii 7.1.3 MeiModel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 7.1.4 Piechowski Model . . . . . . . . . . . . . . . . . . . . . . . . . 105 7.2 Selection and Summary of FinalModel . . . . . . . . . . . . . . . . . 106 8 Model Implementation in EnergyPlus 107 8.1 Buried PipeModule Description . . . . . . . . . . . . . . . . . . . . . 107 8.2 EnergyPlus Code Changes . . . . . . . . . . . . . . . . . . . . . . . . 109 8.2.1 WeatherManager . . . . . . . . . . . . . . . . . . . . . . . . . 109 8.2.2 PlantManagers . . . . . . . . . . . . . . . . . . . . . . . . . . 115 8.2.3 Plant Loop Equipment . . . . . . . . . . . . . . . . . . . . . . 115 8.3 Error Handling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 8.4 Preliminary Evaluation of theModel . . . . . . . . . . . . . . . . . . 116 8.4.1 Input File Description . . . . . . . . . . . . . . . . . . . . . . 116 8.4.2 Discussion of Results . . . . . . . . . . . . . . . . . . . . . . . 117 8.4.3 Direct Implementation Comparison: EnergyPlus vs. VBA . . 121 8.5 EnergyPlus Documentation . . . . . . . . . . . . . . . . . . . . . . . 123 8.6 Computation Time . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 9 Conclusions 125 9.1 Model Development . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 9.2 FinalModel Implementation . . . . . . . . . . . . . . . . . . . . . . . 126 9.3 FutureWork . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 A Earth Tube Model Code 138 B Gambit Journal File 143 C Entering Water Temperature Summary 145 D Modified Piechowski Model Code 150 viii E Modified Mei Model Code 168 F Analytic Constant Wall Temperature Code 185 G EnergyPlus IDD Object 187 H EnergyPlus IDF Object 189 I EnergyPlus Input Output Reference Contribution 191 J EnergyPlus Engineering Reference Contribution 195 ix LIST OF TABLES Table Page 2.1 Variable Definitions for equation 2.1 . . . . . . . . . . . . . . . . . . . 6 3.1 Results of Grid Independence Study . . . . . . . . . . . . . . . . . . . 26 4.1 Earth Tube Model Inputs . . . . . . . . . . . . . . . . . . . . . . . . 30 4.2 Variable Descriptions for Eq. 4.1 . . . . . . . . . . . . . . . . . . . . 31 4.3 Variable Descriptions for Eq. 4.7 . . . . . . . . . . . . . . . . . . . . 32 4.4 Calculation of Earth Tube Outlet Temperature . . . . . . . . . . . . 33 4.5 Common Input Data For Earth Tube Models . . . . . . . . . . . . . . 36 4.6 Convection Correlation Information for Circular Pipe Flow . . . . . . 42 4.7 Earth TubeModel Sensitivity Study Results . . . . . . . . . . . . . . 47 4.8 Earth Tube Model Influence Coefficients and Percent Error . . . . . . 48 5.1 MeiModel Update Equation Organization . . . . . . . . . . . . . . . 54 5.2 MeiModel Nomenclature . . . . . . . . . . . . . . . . . . . . . . . . . 55 6.1 Symbols and Variables used in Piechowski Model . . . . . . . . . . . 80 6.2 Piechowski Model Symbol Definitions . . . . . . . . . . . . . . . . . . 83 6.3 Piechowski Input Variable Tabulated Results . . . . . . . . . . . . . . 87 6.4 Simulation Parameters for Piechowski Verification . . . . . . . . . . . 98 6.5 Piechowski Verification: Exiting Temperature RMS Error . . . . . . . 99 6.6 Piechowski Verification: Heat Transfer Rate RMS Error . . . . . . . . 101 8.1 Computation Time Results . . . . . . . . . . . . . . . . . . . . . . . . 124 x LIST OF FIGURES Figure Page 1.1 Three Possible Horizontal Pipe Configurations . . . . . . . . . . . . . 2 2.1 Kusuda and Achenbach Ground Temperature Prediction . . . . . . . 20 3.1 Generic Fluent Buried PipeModel . . . . . . . . . . . . . . . . . . . 25 3.2 FluentModelMesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 4.1 Earth TubeModels: Exit Air Temperatures . . . . . . . . . . . . . . 38 4.2 Earth TubeModels: Exit Air Temperature Deviations . . . . . . . . . 39 4.3 Earth Tube and Fluent Heat Transfer Rates . . . . . . . . . . . . . . 40 4.4 Earth Tube PipeModel Results . . . . . . . . . . . . . . . . . . . . . 43 4.5 Fluent and Earth Tube PipeModel Exit Temperatures . . . . . . . . 44 4.6 Fluent and Earth Tube PipeModel Delta Exit Temperatures . . . . . 44 5.1 MeiModel Domain Description . . . . . . . . . . . . . . . . . . . . . 51 5.2 MeiModel Finite Difference Grid . . . . . . . . . . . . . . . . . . . . 51 5.3 MeiModel Update Equation Locations . . . . . . . . . . . . . . . . . 53 5.4 MeiModel Description: Angular Grid Locations . . . . . . . . . . . . 59 5.5 MeiModel Test Results: ConstantWall Temperature . . . . . . . . . 65 5.6 MeiModel Test Results: Initial Condition Independence . . . . . . . 67 5.7 Mei and Fluent Model Boundary Differences . . . . . . . . . . . . . . 69 5.8 High Flow Test: Mei and Fluent Cross Section Temperature Profiles . 70 6.1 Large Scale View of Piechowski Domain Cutaway . . . . . . . . . . . 74 6.2 Small Scale View of Piechowski Domain Cutaway . . . . . . . . . . . 74 xi 6.3 Case Layout and Node Descriptions for Piechowski Model . . . . . . 77 6.4 Piechowski Radial Grid Distribution . . . . . . . . . . . . . . . . . . 81 6.5 Piechowski Input Variable Study Results . . . . . . . . . . . . . . . . 88 6.6 ConstantWall TemperatureModel Results . . . . . . . . . . . . . . . 93 6.7 Piechowski Initial Condition Independence Results . . . . . . . . . . . 94 6.8 High Flow Rate Test Results: Piechowski Model Cross Section Temperature Distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 6.9 High Flow Rate Test Results: FluentModel Cross Section Temperature Distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 6.10 Piechowski Verification: ExitingWater Temperature Results . . . . . 100 8.1 Main Pipe SimulationManager Routine Flow . . . . . . . . . . . . . 108 8.2 Pipe Initialization Routine Flow . . . . . . . . . . . . . . . . . . . . . 110 8.3 Main Pipe Calculation Routine Flow . . . . . . . . . . . . . . . . . . 111 8.4 Buried Pipe Calculation Routine Flow . . . . . . . . . . . . . . . . . 112 8.5 Buried Pipe Radial Region Routine Flow . . . . . . . . . . . . . . . . 113 8.6 Buried Pipe Water Boundary Routine Flow . . . . . . . . . . . . . . 114 8.7 Results Section 1: 12Feb to 28Apr . . . . . . . . . . . . . . . . . . . 118 8.8 Results Section 2: 31Mar to 8Apr . . . . . . . . . . . . . . . . . . . 119 8.9 Results Section 3: 31Mar to 28Apr . . . . . . . . . . . . . . . . . . 120 8.10 Results Section 4: 2Jun to 17Jun . . . . . . . . . . . . . . . . . . . 120 8.11 Implementation Comparison, EnergyPlus to VBA: Exiting Temperature121 8.12 Implementation Comparison, EnergyPlus to VBA: Total Heat Transfer Rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 xii CHAPTER 1 Introduction Heat transfer from buried horizontal pipes has been studied for a majority of the twentieth century, encompassing such fields as residential and commercial ground source heat pump system design, industrial pipelines, district heating and cooling, buried power cables, and the use of earth tube heat exchangers. These applications have led to studies of soil properties, effects of moisture on underground heat transfer, effects of multiple buried pipes along with different insulation methods. The results of all these studies have shown that the problem is not trivial. There are many factors which increase the difficulty of finding an accurate solution to the problem of the horizontal pipe. The ground is heterogeneous; soil conductivity varies with temperature, soil type, and moisture content. The ground surface changes constantly based on the annual and daily temperature cycles, ground cover materials (grass, concrete slab, etc.), rain, and snow cover conditions. The fluid in the pipe will not maintain constant behavior along the entire length of the pipe, with properties such as viscosity and density changing with temperature. The geometric complexity is increased when multiple pipes are laid in the same trench. Even within this multiple pipe layout there are different configurations, as shown in figure 1.1. These complexities make finding a quick solution nearly impossible. As with any engineering analysis, the solution begins with a set of valid assumptions. The inclusion of a buried pipe model in an annual simulation program is a necessity for properly modeling horizontal ground loop heat exchangers. It is also important for district heating and cooling systems. With lowenergy houses and buildings using 1 Figure 1.1: Three Possible Horizontal Pipe Configurations ground heat transfer, the model becomes even more important, as every degree of precooling required of a water source and every degree of heat lost from a heating pipe can be costly. The main purpose of this current study is to identify a buried pipe heat transfer model which balances accuracy and computation time, and implement that model in the EnergyPlus (2007a) program. The study accomplishes this by the following procedure: 1. Examine previous work in the field. This includes experimental work, which would help define the important aspects of a buried pipe system, and modeling work, which would help define the assumptions that other studies have made in order to form a tractable solution. 2. Develop a reference model. Without experimental data, a robust model must be utilized to provide a comparison tool. 3. Evaluate currently implemented models in EnergyPlus. This step includes looking into models which may not initially be included as buried pipe models, but may be modified to fit such a problem. 4. Develop models for analysis. Each model that seems like a useful option is developed in some form and tested to determine the quality of the results. 2 5. Select a single buried pipe model to implement. The best choice for the model will balance accuracy and computation time, while providing a general purpose solution to fit all applications of buried pipes. 6. Implement the model in EnergyPlus. The model will act as its own module in the EnergyPlus environment, and it will be tested thoroughly to ensure proper operation. 7. Provide user documentation for the model. This documentation will be included as part of the EnergyPlus documentation state what inputs and outputs are available for the model, as well as a reference showing how the model works. During many of the tests and simulations included in this work, hourly weather and entering water temperature data are used as boundary conditions. The weather data are annual hourly data from Atlanta, GA. The entering water temperature data are also hourly data, but it is from an experimental source, the ground source heat pump plant at Oklahoma State. Even though this data is not applicable as model validation, it does provide realistic inputs to the models. This report includes a total of 9 chapters, including this introduction as chapter one. Chapter two is a literature review of different modeling techniques and other aspects of buried pipe design. Chapter three is an introduction to the Fluent reference model used to evaluate the current modeling work. Chapter four evaluates an existing EnergyPlus model, the earth tube model, for modification and use as a buried pipe model. Chapter five evaluates a fully radial model with variable boundary conditions. Chapter six evaluates a model which mixes both cartesian and radial coordinates within the same domain. Chapter seven summarizes results from each model, and selects a model for EnergyPlus implementation. Chapter eight discusses implementation of this model in the EnergyPlus program. Chapter nine summarizes the conclusions of the study and proposes future work. 3 CHAPTER 2 Literature Review In order to find a proper model to implement, it is necessary to look into literature to see what attempts have previously been made. When looking for a buried pipe model, the search may be broadened slightly as well, to include buried ’earth tubes’, since a detailed analysis of a buried air tube shares many aspects with a buried water pipe. The key when looking at these models is to understand the underlying assumptions and ensure that they remain valid for the system at hand. There are several approaches for solving the buried pipe problem, including semianalytic solutions, finite difference solutions, and finite volume solutions. All of these calculation methods are important, but an understanding of experimental work is required in order to validate the applicability of any model. The current search is broad enough to include the experimental studies and understand the important aspects of a buried pipe system design. Some of the other related material comes from district heating literature, which typically includes substantially long piping systems; soil property prediction, which is an important boundary condition for any simulation dealing with the soil; and transport delay, which can be a significant factor, especially with longer piping systems. 2.1 SemiAnalytic Solutions The semianalytic solutions typically consist of a line source or cylinder source solution. In this solution, the soil is assumed infinite in all directions, and an infinitesimally small heat source is applied which represents the pipe. In some cases, this 4 heat source is the input to the model, resulting in pipe outlet temperatures or ground temperatures. The model can be reformulated to instead determine the heat rate when temperatures are given as inputs. Even with this formulation, there are other limitations of the model. The model does not account for ground surface interaction, does not readily account for initial ground temperature gradients, and the line source is also variable along the length of the pipe due to the fluid temperature changing, so this method cannot suitably predict the behavior along the pipe. In many cases these models also use annual average data which can produce a typical annual result. One of the earliest line source models was introduced in a 1948 paper in ASHVE by Ingersoll & Plass (1948). The key in this paper is describing the temperature distribution around a buried pipe. An equation is offered, (2.1), based on the assumptions mentioned above, that can give the temperature of any radial distance away from the line source, valid from the pipe wall to a large distance away. This model assumes that the pipe is the only cause of temperature change in the soil. T − TO = Q 2πk ∞ x e−β2 β dβ = Q 2πk I(x) (2.1) where the variables are defined in table 2.1. The paper offers values of I(x) for numerous values of x, and an equation approximating I(x), which would be more suitable to a simulation program than a lookup table. The original equation was only valid for a small range of x. For the purpose of this current thesis work, the equations were adjusted, and equation (2.2) results, and is valid for a wider range of x. I (x) = 2.303 · log 1 x + x2 2 − x4 8 + 0.1128 · x3 − 0.2886 (2.2) The rest of the paper deals with different pipe situations, and solving examples. The examples range from simple single buried pipes to double buried pipes, and then to simulations where the heat rate changes from month to month. Although these 5 Table 2.1: Variable Definitions for equation 2.1 Variable Definition x √r 2αt T Temperature in soil at any selected distance from the pipe [deg F] To Initial temperature of soil [deg F] Q’ Heat emission of pipe (negative for absorption) [Btu/linear ft/hr] r Distance from centerline of pipe [ft] k Thermal conductivity of soil [Btu/hr/ft2/F] α Thermal diffusivity of soil = k ρCp ρ Density of soil [lb/ft3] Cp Specific heat of soil [Btu/lb/F] t Time since start of operation [hr] β Variable of integration examples do show that the line source approach is a widely applicable method, it still fails in the limitations discussed previously (uniform boundaries, no ground surface interaction). The method still offers incredible simplicity, and has been used for a significant time frame, but with simulations requiring more accuracy, especially in green design, this method does not seem to be the ideal solution. The simplifications for line source applications can be modified and built upon to account for different pipe configurations. In Persson & Claesson (2005), the multipole method is combined with the line source approach in order to account for multiple pipes using a steady state analysis. In this situation, each pipe is considered as a line source, and the resulting temperature distribution around the pipe can be described by an infinite series solution. When this is performed for each pipe, the temperatures and heat flows can be constrained to ’line up’ in between the pipes, thus creating an infinite series of equations describing the system. This infinite series can be truncated 6 to a desired level of accuracy, resulting in a solution accounting for both pipes. In this paper, multiple pipes are simulated in a vertical configuration (one pipe above the other), as well as a horizontal configuration (each pipe next to the other). In the vertical configuration, the supply pipe is simulated both on top and on the bottom. The pipe system is modeled using district heating data (geometry, temperatures) as inputs. It is found that using these simplifications, a fourth order approximation to the infinite solution can provide results of significant accuracy with a computation time of a few seconds. A different analysis using line sources is found in Saastamoinen (2007). In this analysis, a rigorous development of unsteady state equations is performed. The application is a floor heating system, in which there is a conducting layer within the floor where heating pipes are contained. Above this, a layer of carpet is modeled, and underneath the conducting layer is insulation. The development includes a convective boundary condition on the top, and constant temperature boundary condition on the bottom. The unsteady behavior is modeled as a fourier series of periodic onoff cycles. The model assumes the shapes of the isotherms around the line sources are nearly circular near the pipe in order to simplify the model to a tractable solution. The effects of ranging placement of the coils and surface Biot numbers are studied in steady state mode and discussed. Another approach which uses semianalytic methods is performed by Arimilli & Parang (1983) on dual horizontal buried pipes. The method used is a boundary integral method. In this method, a steady state version of Laplace’s equation, represented in Eq (2.3), is written for a twodimensional domain according to boundary integral methods, with the modified version given in Eq (2.4). ∇2T = 0 (2.3) 7 1 2 T (Y ) + B T ∂g ∂n − g ∂T ∂n ds = 0 (2.4) In Eq (2.4), the B subscript refers to integration over the entire boundary, T(Y) represents the temperature at any point Y on the boundary, and g(X,Y) is a harmonic function over the solution domain. The boundary integral involves the temperature along the boundary as well as the derivative at the surface (heat flux for the thermal analysis case). By using this method, the integral can be discretized over a set of segments in the domain. This is similar to a finite difference method, however the boundary integral approach only involves values at the boundaries, and not intermediate values in the center of the domain, which are not of interest in most studies. The resulting equations can be solved with iterative means until convergence. This paper provides the equations required to solve for the general case and describes their simplicity of use, although no validation or verification is given so the accuracy is uncertain. 2.2 Numerical Analysis In many cases, the geometry may be too complex to model properly with a line source approach. For these cases a numerical analysis lends itself well. The numerical approach could be a simple 1D finite difference model, or a 2D finite difference model with variable boundary conditions, or a full 3D finite volume model with unstructured boundary conditions. The approach taken would depend on the pipe configuration, accuracy required, and computation time allowed. The following works span a broad spectrum of numerical solutions and each can provide their own insight into a satisfactory model for implementation in an annual simulation program such as EnergyPlus. There are many variations in ground heat exchangers including concentric tubes, horizontal trenches, and vertical boreholes. In order to develop a model of a vertical heat exchanger, some of the same assumptions can be made, and some different 8 assumptions also. For instance, in a vertical heat exchanger, there will almost always be two tubes in very close proximity. In this configuration, it would be a very bad approach to neglect interaction between the pipes. In many cases, a single horizontal heat exchanger will be placed in a trench, in which case this assumption would be valid. For a vertical heat exchanger, the effects of ground surface will be minimal, as the majority of the pipe will be buried deep underground, while in a horizontal heat exchanger, the entire length of the pipe is affected by the behavior above the ground surface. Chiasson & Spitler (2000) describes the model development of a heated bridge deck using a ground source heat pump system. Lei (1993) presents a validated, computationally efficient model for vertical ground heat exchangers. The model is simplified by neglecting the interface resistance between the pipes and the surrounding material. The effects of the Utube at the end of the pipe are also neglected. The model’s computational efficiency is improved by using a dual domain configuration. A fine region is used in the high gradient regions near the pipe, along with a fully implicit solution scheme. The coarse region farther away is solved using an explicit scheme. The time step used is efficiently selected based on the geometry. In both regions, adaptive grid spacing is employed to avoid wasting computation power where it is not needed. With these assumptions and improvements in place, the cylindrical transient equations are discretized and applied to each node in the domain, and the simulation is solved. The results are compared to experimental results, with high accuracy while the piping system is active. When the flow is stopped, as in the pumping off cycle, the results deviate somewhat. The exiting water temperature during these time steps is not drastically important in most engineering analysis, as long as the model simulates correctly once the flow comes back on. A set of studies is performed by Mei (1988, 1986a) in which both single and double horizontal pipes are simulated. In Mei (1988), a dual pipe is studied, with particular interest in the effects that the second pipe has in the temperature distribution around 9 the first pipe. A quasi3D finite difference code is developed, based on energy balance equations, which simulates both a single pipe and a dual pipe in a tubeovertube configuration. The model assumes that the temperature from inside the pipe up to the pipe wall is radially symmetric, and a circular far field boundary condition is applied which is centered in the bottom tube, and with the top at the ground surface. In order to simplify the model, axial heat transfer is neglected in the pipe and soil, and therefore axial heat is only allowed to flow through the fluid. The model results are compared to experimental data in terms of daily energy absorption in the ground, with the model producing an average error of less than 12%. Mei (1986a) describes three models, each with their own distinct purpose: soil freezing, seasonal ground temperature variations, and thermal interference. Each model is described in detail including mathematical development, computer implementation, validation, and parametric studies. It is found that for soil freezing around the pipe, the effect is minimal as long as the entering temperature is not less than 4◦C, which is typically used as a design constraint anyway to maintain sufficient heat pump capacities. The seasonal ground temperature variation model provides a very high quality simulation of the temperature distribution around a pipe, and also allows for studies of different backfilling materials. This model seems promising for simulation applications and is evaluated further in chapter 5. The thermal interference model is used for any case with multiple buried pipes in the same trench, and is not investigated in the current work, although analysis of multiple buried pipes is recommended as future work. A study which involves both heat and mass transfer was performed by Piechowski (1996, 1998, 1999). In the thesis, Piechowski (1996) described every aspect of the model development in great detail, while the other two papers, Piechowski (1998, 1999), give insight into the model validation and theoretical development, respectively. The model contains aspects that most buried pipe models do not, such as moisture diffusion through the soil, and the addition of thermal storage into the system. The 10 model assumes that the soil temperature a certain distance away from the pipe is not affected by the pipe presence. This becomes the far field boundary condition, and is based on a diurnal and seasonal temperature variation. In order to model moisture transfer in a tractable manner, the mass transfer is only modeled in the vicinity near the pipe, which is where the highest gradients will exist for both heat and mass transfer. The moisture equations therefore have a boundary condition which is in this vicinity. The heat transfer is simulated over a larger domain, with variable boundary conditions. The system is simulated in a 3D finite difference domain, and as with other models, heat and mass transfer are neglected in the axial direction in the pipe material and the soil. This allows for a set of 2D domains to be solved, connected only by the fluid transferred through the pipe. The model is validated against experimental data, and it is found that for a small temperature range, the inclusion of moisture effects (latent heat and moisture migration) does not significantly increase the accuracy, so for most models is would not be necessary. However, the soil thermal conductivity is a very significant input, and improper assumptions of soil properties will lead to inaccurate results. The Piechowski model is developed and tested in chapter 6. Tarnawski & Leong (1993) presented a model of a ground coupled heat exchanger, but with more emphasis on the overall system simulation investigated. The model includes data for several heat pumps, soil properties, heat exchanger configuration (spiral HX included), and pipe materials. The model is verified against simulated results and previously published experimental work, and agrees with a fair amount of accuracy. The verification process was hindered by a lack of information from the verification source, including ground properties and initial conditions. 2.3 Response Factor Models Another heat transfer model type involves the use of transfer functions or response factors. Transfer functions are typically confined to certain systems while response 11 factors can be developed for any system. These are already found in building simulation programs as a way to calculate transient heat transfer through walls using pseudo steadystate means. For each wall construction, transfer functions or response factors are developed which represent how the wall responds if there is a one degree change in temperature on one surface. Then superposition can be applied to determine the effects of general temperature changes. There are also response factors used in ground heat transfer, but there are limitations to this method. A short time step response factor model for vertical borehole heat exchangers was developed by Yavuzturk & Spitler (1999) and validated by Yavuzturk & Spitler (2001). In this model, a series of response factors are developed using a transient two dimensional finite difference system. With the response factors calculated, the model is able to produce the temperature response of the borehole field to a step heat input from the ground heat exchanger. This model type is appealing because once the response factors have been generated, the responses are computationally fast. The response factors developed in this work are also valid for short time steps, which is suitable for an energy simulation program such as EnergyPlus. This model proves to work quite well for vertical heat exchangers. Horizontal heat exchangers have a disadvantage for such an approach in that they are greatly affected by the ground surface. For a vertical heat exchanger, only the small section at the top will be affected, so the response factors are allowed to ignore this behavior. For a shallow horizontal heat exchanger, the response factors would need to include the response in the ground to not only a step heat input from the heat exchanger, but also to a step change in temperature or heat from the ground surface. The response factors must also be able to handle modeling a certain depth of soil. Currently, there is a horizontal surface heat exchanger model available in the EnergyPlus program that is based on a modified model from Strand (1995). This model uses transfer functions which are modified to include a heat source within the construction. This is also found when modeling radiant floor systems, in which 12 case a warm pipe is embedded within the construction. This method is quite useful and again, it is very fast compared to iterative finite difference techniques (the transfer functions in place are also iterative due to required surface heat balances). The major problem with using this model as a general purpose buried pipe model is the depth. Existing transfer functions do not work with very heavy constructions, and with approximately 1m of soil, the EnergyPlus transfer functions diverge. Response factors can be developed for any system such as a horizontal heat exchanger, however the work is beyond the scope of this thesis. The development of response factors is recommended as future work. 2.4 Experimental Work To ensure that any model is relevant and useful, there must be some experimental work in the field. This experimental work can provide a means of validating a model, provide ranges of useful input data, characterize important properties of the system that must be investigated and accounted for in the model, and also provide a proof of concept to ensure that model development is worth the time and effort. Inalli & Esen (2004) provided details on a horizontal ground source heat exchanger system in heating season in Turkey. The study provides key data which could be used as input to a new model development. The study also provides experimental evidence that the soil properties play a significant role in the design of the heat exchanger, as it can easily lead to oversizing of the components and therefore the whole system. In this study, the average coefficients of performance (COPs) of the system (2.66 2.81) were lower than in other experiments due to oversizing of the equipment. The data provided is presented more in terms of seasonal and daily averages. The data provided is not suitable as a source for model validation. An article that discusses not only experimental apparatus, but also useful design techniques, was presented by Klimkowski et al. (1985). The paper studies the long 13 term effects of average Uvalue for a given horizontal heat exchanger. Data is measured, and then resolved back into average conductance. After a period of continuous operation of the heat pump, the Uvalue approaches a steady constant, one value in winter, and another value in summer. It is recommended to use this value as a design criteria. The paper then offers a design process, which simply utilizes average earth temperature, fluid temperature, and a given Uvalue, along with a required capacity, to calculate a required HX length in both summer and winter. The longer of these two results is selected as the required pipe length. This type of analysis is a very simplified model which may provide a decent prediction of total annual energy usage, but gives no insight into the daily or hourly operation. These smaller time steps, along with accounting for energy storage in the ground, is where savings may be found using a more robust model, especially when looking at low energy applications. 2.5 Earth Air Heat Exchanger Analysis When developing a model for a buried pipe, it is desirable to look to different aspects of the heat transfer field in order to possibly find a suitable model that can be modified to fit the needs of a buried pipe. One example is an earth tube. In an earth tube, air is passed through a buried duct in order to precool the air in summer, and preheat the air in winter. In this way, the demands of the conditioning equipment are reduced, while the ventilation requirements are easily met with fresh air. The earth tube models differ in the assumptions, as air is passing through the pipe instead of water or glycol. However, if a model shows good performance under earth tube conditions, it may be useful to modify it and analyze it further as a buried pipe. AlAjmi et al. (2006) presented an earth tube model which employs a variable soil temperature variation, and applies steady state behavior and an effectiveness correlation in order to calculate exiting temperature conditions and overall heat transfer. The model is verified against previously published data and proves a high level of 14 accuracy. The model is then implemented into the TRNSYS environment and used in a full building simulation to show the benefits of using an earth tube, as the cooling load is decreased substantially. Baxter (1992, 1994) provided a full description of experimental results of an earth tube system in both heating and cooling modes. The results show that one of the main factors that affects the earth tube performance is the temperature gradient that forms in the axial direction near the pipe. This indicates that assuming constant soil temperature near the pipe should be used cautiously, as it could be erroneous. Bojic et al. (1999) simulated a pvc and a steel pipe together, and showed results of the model. The model was based on a 3D Finite Element Method, although no details were given about the meshing specifics. While the results are given, no validation or verification is offered. One conclusion is that the earth tube operates better in the summer months than the winter months. This may be due to the simulation climate and other temperatures used in the simulation. Another approach to modeling an earth tube is given by De Paepe & Janssens (2003). In this paper, the authors developed a thermo hydraulic design of earthair heat exchanger. In this model, both the heat transfer effects and the fluid mechanics effects are modeled, resulting in a very detailed study. The purpose of this model is to evaluate the balance between pressure drop and heat transfer when designing a heat exchanger. The paper develops an NTU correlation for the earth tube heat exchanger. A ’specific pressure drop’ is developed to allow a designer to then calculate required lengths of parallel tubes in order to achieve the required capacities. A very simplified model of an earth tube, which is implemented in the EnergyPlus simulation program, is developed by Lee & Strand (2006). This model is developed in a very simplified manner in order to achieve two main goals. The first goal is implementation in an annual simulation program. A full finite volume study of an earth tube would be too computationally intensive to recalculate every time step 15 of an annual simulation. The second goal is available input data. A simulation engine such as EnergyPlus requires that the models have input to which designers have access. The model takes in simple average annual data, and at each time step, calculates a soil temperature in the near pipe region, and calculates soil, pipe, and convective resistances. A steady state approach then uses an energy balance between the soil and air to create an equation to calculate exiting air temperature as a function of pipe geometry, fluid characteristics, and soil data. The model can also account for fan power if mechanical ventilation is selected. This model is an almost instantaneous calculation model, and therefore lends itself well to being converted to a buried pipe model and reanalyzed. This model is evaluated in detail in chapter 4. Mihalakakou et al. (1994) provides a model which discusses both heat and mass transfer around the pipe, and accurately predicts the exiting air temperature. The model uses a 3D finite element approach which uses concentric rings as the elements. The model offers the conservation equations which are applied to each node and include both heat transfer and moisture transfer. With boundary conditions specified, the model was used to simulate an earth tube over a period of 15 days. Over this time, experimental data was measured, and the results of the model were in very satisfactory agreement. This model was implemented within the TRNSYS environment, and provides a means of properly predicting earth tube exiting temperature as well as soil temperature variations. In a study by Tzaferis et al. (1992), two sets of experimental data were measured, and eight different earth air heat exchanger models were simulated under these experimental conditions, and their results were compared, along with the model sensitivity. The models were split in two categories: those which analyzed ground heat transfer, and those which only analyzed heat transfer through the pipe. In the latter of these, temperatures at or near the pipe wall were required as inputs. The models range from assuming constant pipe wall temperature (analytic integral) to using finite difference 16 schemes. The root mean square (RMS) error for nearly all the models were below 3.5%, showing that even simplified models can provide decent results for an earth tube. 2.6 District Heating One of the larger scale applications of a buried pipe model is in district heating systems. In these systems, a long pipe loop connects a large number of buildings, providing hot water, which can then be tapped and used as a primary heat source or as a source for heat pumps. In these systems, a good model for the pipe heat transfer effects should be used since even a small percentage of energy savings can be quite significant. A simplified model by Bohm (2000) takes an approach quite different from other models. The model assumes that at each time step, there is a soil temperature that can be used in steady state equations to mimic transient behavior. This temperature would account for the storage effects in the ground by using a temperature from a previous time step. In this manner, the soil temperature near the pipe would be ’lagged’ by a certain amount. The position of this undisturbed temperature in the ground is found by numerical simulations and experimentation. A function is developed for this temperature, although several restrictions were placed on the analysis. For the verification, a constant pipe temperature was used, and it is recommended to only maintain constant model inputs for very brief periods. In this study, as with most district heating systems, the pipe configuration is a dual pipe, with the pipes laid side by side. The results do compare well with experimental data taken, although, some of the model inputs come from predicted, or experimentally determined data. A model of an indirect district heating system is offered by Chuanshan (1997). This model looks deeply at very detailed effects in a district heating system, and not so deeply at predicting the heat transfer of the buried piping. The study does offer 17 details on model parameters which are very useful if a buried pipe model is developed and tested in a district heating simulation. Eriksson & Sunden (1998) detailed some new aspects of modeling, while still accounting for both heat and mass transfer in the soil. The model appears to be a 1D approach, with the solution using a tridiagonal matrix algorithm (TDMA) approach. The key to this particular model is the investigation of pipe material breaking down through time. The thickness and conductivity of the pipe changing over time is studied over a longterm simulation, and the results show that the thermal conductivity of the pipe increases over time, but does not make a significant impact in the overall heat exchanger performance. 2.7 Soil Property Prediction In order to properly simulate a ground coupled heat exchanger, the soil effects must be modeled accurately as well. The effects of surface temperature variation is even more important when dealing with a horizontal heat exchanger than with a vertical heat exchanger because the pipe or duct is near the surface for the full length, not just a small initial section. This section discusses models and information found for predicting not only ground temperatures, but also ground properties. A study by AbuHamdeh & Reeder (2000) shows that the thermal conductivity is highly dependent on other soil properties such as density, but the application of the study is not relevant for a buried pipe calculation. Fischer (1983) gives a brief summary of fourteen soil heat and moisture transfer models, describing each model’s capabilities and limitations, but does not include any model development itself. A model in Kusuda & Achenbach (1965) is one of the more popular models for predicting ground temperatures. The model is based on an assumption that the earth is a semiinfinite solid, and that the surface boundary condition can be described by a sinusoidal temperature variation. In this manner, a simple correlation is given, based 18 on some simple input data, but the ground surface temperature and average annual amplitude of temperature variation must be known. The main equation (2.5) uses an exponential term to account for depth into the ground (delay), and a cosine term to account for annual variation. T = A − BO · e −x· √ π DT cos 2πθ T − x π DT − PO (2.5) In equation (2.5), T is the ground temperature at a given time θ and depth x. D is the thermal diffusivity of the soil, T is the time period of the annual variation, and should be the same units as θ. PO would be a phase angle to account for the date of minimum surface temperature. A is the annual average surface temperature, and BO is the annual average variation in surface temperature. As an example, to show the effects of ground heat diffusion, figure 2.1 shows the annual ground temperature for a series of depths. In this example, A was set to 18.333◦C, BO was set to 11.667◦C, D was set to 1.765E3 m2/hr, PO was set to zero, and since the run was an annual simulation and diffusivity was in hours, the constant T was set to 8760 hours. Kumar & Kaushik (1995) discusses a soil temperature model which is based on Kusuda & Achenbach (1965), but includes soil properties (conductivity, specific heat, density, and diffusivity) for ten different soil consistencies (dry light, damp heavy, average rock, etc.). 2.8 Transport Delay One aspect that is typically not dealt with in detail is the transport delay problem. Most simulation programs have pipe models, but do not account for the time required for the fluid to pass through the pipe. In the past, this may have been neglected due to the additional computational strain that this would place on a simulation program. With computers becoming faster, this aspect should not be ignored when it is applicable to include it. 19 Figure 2.1: Kusuda and Achenbach Ground Temperature Prediction 20 Hanby et al. (2002) develops a model that is based on second order equations used at several nodes along the length of a pipe. The model accounts for the thermal capacitance of the fluid and the conduit wall. Each node is a wellmixed zone, and the node separation distance is based on capacity values and average fluid residence time. This residence time is a factor which can be assumed based on flow behavior or precalculated using CFD studies. The model shows very close agreement with experimental data. Saman & Mahdi (1996) details some previous models that attempted to predict the transport delay problem, along with where these models fell short. The model in the paper is based on analytic energy balance equations discretized along the length of the pipe, and solved using finite difference methods. The results do follow the expected behavior, but no validation is offered in the paper, so the accuracy of the model cannot be determined. 2.9 Literature Summary This chapter has been a review of literature involving buried pipe heat transfer models, soil property prediction, and transport delay models. The buried pipe models included semianalytic methods, numerical simulation methods, and response factor methods. The response factor methods reviewed were generated for vertical heat exchangers, which allowed the ground surface to be ignored. The generation of response factors for a shallow horizontal heat exchanger is not in the scope of this thesis work, but is highly recommended as future work. The semianalytic methods required a large number of assumptions and typically required input such as the line source strength; doing so requires a guess and therefore introduces high uncertainty into the results. The numerical methods varied in computational requirements and accuracy. The models selected for development in this work are the Earth Tube model, the V. C. Mei model, and the Piechowski model. These three 21 give varying degrees of assumptions and accuracy, and should give a nice comparison of what work has been performed, and which model is worthy of development within the EnergyPlus system. The soil property prediction by Kusuda has been used for years covering different fields of study. For the models currently developed in this work, the Kusuda and Achenbach method will be used to predict boundary condition temperatures and surface temperatures when appropriate. The transport delay models will currently not be developed, as the main goal of this thesis work is to get a buried pipe heat transfer model into EnergyPlus. The introduction of a transport delay model could be obtained by linking several smaller pipe models together once one is developed. This is left as another good recommendation for future work. However, with a transient buried pipe model, the effects of a transport delay are somewhat captured. The basic pipe model allows the inlet temperature to pass directly through to the pipe outlet. With the addition of any transient effects, the immediate transport will encounter a delay. Further investigation will reveal if the models need to be developed together to create one dualpurpose model. In conclusion, three numerical heat transfer models will be developed which use the Kusuda and Achenbach soil temperature prediction method, and do not include an actual transport delay aspect. These three models will be compared against a verification model discussed in chapter 3, and a decision will be made to implement one of the models weighing accuracy vs. computation time. 22 CHAPTER 3 Fluent Reference Model 3.1 Fluent Model Description When creating or testing any model, there must be a method of ensuring the model is operating properly. In some cases, analytic solutions may be available which can be used as a validation source. Another way of performing validation is through experimental results. If one can create an experimental apparatus or facility to which the model can be applied, the results can be compared to prove the model’s validity. In the current case, experimental data is not available, and so direct validation is not possible. In some limiting cases, analytic means are used in this work to prove model operation. These cases require substantial assumptions to allow an analytic solution to be developed. For the majority of this current work, verification will be performed. In verification, the model is compared against a more robust, or previously validated model. Fluent CFD software is employed for the current situation. In most cases, Fluent is used to simulate fluid flow in diverse situations. However, the finite volume engine within Fluent is powerful enough to simply switch off the momentum equations, and only solve the conservation of energy equation in solid regions. This allows conduction to be modeled easily within Fluent, along with convection flow in a tube, all in the same model. Another benefit of using Fluent is the flexibility of boundary conditions. There are features within Fluent that allow users to input their own functions that can be spatially or temporally dependent, and applied to any boundary of the model. The functions are written in the C programming language, can be easily hooked to a 23 model, and can access many internal variables to Fluent’s solver. This is very useful when simulating, for instance, the annual variation of ground surface temperature. The user defined function can reach within Fluent, retrieve the current simulation time step, and return a value of surface temperature for that particular hour of the year. For some situations where the boundary condition is spatially dependent, such as the Kusuda and Achenbach boundary (a function of time and depth), the function can reach into Fluent and get the depth coordinate in order to return the appropriate earth temperature for that time and depth. The other boundary conditions can all be applied as needed, some symmetric, some known boundary temperature. The exact configuration will depend on the model being verified. Figure 3.1 shows a typical Fluent model for this work. The buried pipe channel along with the relevant boundary conditions, are shown. The boundary conditions will vary depending on the current model being tested. Throughout this work, a number of simplified models will be compared to Fluent models that have been carefully generated to closely match the simplified model. The details of these comparisons are discussed in each individual model chapter. 3.2 Fluent Model Grid Refinement In order to develop a valid Fluent model, an independent grid must be demonstrated. By using the gambit journal file, it is easy to reproduce models quickly. An example of the journal files used in this study can be found in appendix B. In order to ensure that the grid is independent, the grid was reproduced several times, each with increasing refinement of the grid. The grid was refined in both the pipe and soil regions, each with its own given node spacing value. For the grid independence check, a steady state situation was utilized, with simple boundary conditions. The ground surface was held at 310K, while the other spatial boundaries were adiabatic. The pipe inlet was air at a uniform 1 m/s and 260K. The entire domain was initialized to the conditions 24 Figure 3.1: Generic Fluent Buried Pipe Model 25 Table 3.1: Results of Grid Independence Study at the pipe inlet, and the simulation ran to convergence. Convergence was kept at the default residual values (1E03 for most variables, 1E06 for energy equation). After convergence was attained and the results were checked for satisfactory convergence (no oscillation, etc.), the total heat transfer rate on the pipe wall was calculated by Fluent. This was used as the grid independence monitor. The grid was simulated for four different mesh conditions, and table 3.1 shows the results of the runs. Table 3.1 shows that the results are in fact converging on a solution as the mesh becomes more and more refined. The table also shows that the time required for the simulation to run is becoming much larger. In order to provide quality results in a reasonable time frame, a compromise must be made between accuracy and computation time. For further computations, the mesh for case 3 (10 minute run time) is selected. This simulation gives quite good results with only about 5% improvement over the previous run. The mesh is shown in figure 3.2. The grid independence is one of the most important checks to make when performing finite difference studies. Other studies to verify the Fluent model under different circumstances will be performed on a modelbymodel basis throughout this work. 26 Figure 3.2: Fluent Model Mesh 27 CHAPTER 4 Earth Tube Model 4.1 Earth Tube Model Description One of the models which is used to predict ground heat transfer from a buried pipe was developed for an Earth Tube system (Lee & Strand 2006). In this earth tube system, ventilation air is brought into an underground duct from a certain distance away from the zone or conditioning system, and the air is tempered, either preheated or precooled depending on the season, and then brought into the zone or through the conditioning equipment. This pretreatment of the air helps to reduce the overall heating and cooling equipment loads while satisfying fresh air requirements. The model makes many assumptions, but requires minimal input data and has a very short computation time. These features lend themselves well to first stage design work in an annual simulation program. Because of the many assumptions, this model seems to be one of the simplest models to apply to a buried pipe. As with many models, one must begin by understanding the underlying assumptions, in order to ensure that the model is valid for a particular system, and not oversimplified. • The first assumption is that fluid flow within the pipe is thermally and hydrodynamically developed. This assumption can be easily made when the length of the pipe is much greater than the diameter. Developing flow will only occur within a short initial section of the pipe, and will not greatly influence the thermal behavior when the pipe is long. This assumption allows a single correlation 28 to be used for any given fluid and flow regime. • The next assumption is that soil temperature in the pipe vicinity can be calculated using the model developed by Kusuda & Achenbach (1965). This earth tube model makes use of the soil temperature model to calculate a soil temperature one annulus from the center of the pipe. • The temperature in the soil near the pipe is not affected by the pipe itself. This assumption allows the soil temperature to be uniform along the entire axial length of the pipe. This allows the solution of the system to be greatly simplified, but may impact the accuracy. • The soil is homogeneous and maintains constant thermal properties (thermal conductivity) throughout any time step, although at each time step, the properties can be adjusted to a new temperature dependent value and held constant yet again. • The fluid is axisymmetric at any cross section, and the pipe maintains a constant cross sectional area in the axial direction. • At any given time step, the heat transfer behavior can be modeled as steady state. In order to be useful, the model inputs must be available. Typically, this requires simplifying the model to accommodate data that most people can access. This particular model requires input of properties such as system geometry, material properties, fluid flow properties (including a convection correlation), and annual surface temperature data. The full set of input data is described in table 4.1. The system geometry should be available to a designer, and the fluid flow properties can be approximated based on design conditions. Soil properties and annual surface temperature data must be estimated. The soil properties will be variable based on moisture content, 29 Table 4.1: Earth Tube Model Inputs Pipe Data Geometry: Buried depth, Axial length, Radius, Wall thickness Properties: Pipe wall thermal conductivity Soil Data Properties: Thermal diffusivity, Thermal conductivity Other: Annual average surface temperature Annual surface temperature amplitude Fluid Data Flow: Flow rate, Convection correlation Properties: Density, Specific heat, Prandtl number Viscosity, Thermal conductivity and the surface temperature data will be location specific. Soil properties may be determined from an onsite thermal conductivity test, and surface temperature data can be approximated from annual weather data. Once the model has the required input data, it uses hourly outdoor air temperature and a series of equations to calculate the air temperature exiting the earth tube. These equations are described here. The first equation calculates the ground temperature in the vicinity of the buried pipe. This equation is developed from the work by Kusuda & Achenbach (1965), but arranged to fit the attributes of this particular model: Tgnd,z = Ts − As ∗ e −z∗ √ π 365∗α ∗ cos 2 ∗ π 365 ∗ t − tshift − z 2 ∗ 365 πα (4.1) The variables used in equation (4.1) are defined in table 4.2. 30 Table 4.2: Variable Descriptions for Eq. 4.1 Symbol Description Tgnd,z Ground temperature at pipe depth [C] Ts Annual average surface temperature [C] As Annual surface temperature amplitude [C] z Pipe burial depth [m] α Ground thermal diffusivity [m2/day] t Simulation Run Time [day] tshift Day of minimum surface temperature (i.e. 0 = 1Jan) [day] With the ground boundary temperature found at the pipe depth, focus can be switched to the fluid. To define the fluid to pipe interface, a convection coefficient is required. The correlation (not listed here) is a function of the working fluid (air, water) and the flow regime (laminar, turbulent). For typical fluids such as air and water, there are many correlations available in common heat transfer and fluid flow textbooks. Whichever correlation is used, the outcome is the convection coefficient, h, which can then be used in further calculations. The soil data, pipe data, and fluid data have now all been initialized and steps can be taken to compute the exiting fluid temperature. First, the thermal resistances of each material and interface are calculated, then the total resistance, and therefore the total conductance, can be calculated with the following equations. Convection resistance : Rc = 1 2πr1h (4.2) Pipe wall resistance : Rp = 1 2πkp ∗ ln r1 + r2 r1 (4.3) Soil resistance : Rs = 1 2πks ∗ ln r1 + r2 + r3 r1 + r2 (4.4) Total resistance : Rt = Rc + Rp + Rs (4.5) Total conductance : Ut = 1 Rt (4.6) 31 The paper then applies a heat balance to a differential length of pipe, where the heat transfer to the air is equated to the heat storage within the air. Because of the simplifications that previously resulted in a uniform temperature along the length of the pipe, the outlet air temperature is a simple expression which is a function of conductance, fluid properties, entering fluid temperature, ground temperature near the pipe, and pipe length. For simplification, the paper defines an intermediate term that can be calculated, which simplifies the final exiting temperature expression, and is given as equation (4.7). A = m˙fCf ∗ ln Tf,i − Tgnd,z − Ut ∗ L m˙fCf (4.7) The variables used in equation (4.7) are described in table 4.3. Table 4.3: Variable Descriptions for Eq. 4.7 Symbol Description m˙f Mass flow rate of fluid Cf Specific heat capacity of fluid Tf,i Entering fluid temperature Tgnd,z Ground temperature near pipe L Pipe length Using this intermediate term from equation (4.7), the outlet temperature can be calculated based on different scenarios regarding entering fluid temperature, Tf,i, and ground temperature, Tgnd,z, as described in table 4.4. As this description shows, this model is a very simple model to implement, and lends itself well to annual simulation programs, where it may be called thousands 32 Table 4.4: Calculation of Earth Tube Outlet Temperature Situation Equation Tf,i > Tz,t Tf,ex (L) = Tz,t + eA Tf,i < Tz,t Tf,ex (L) = Tz,t − eA Tf,i = Tz,t Tf,ex (L) = Tz,t of times during one simulation. When implementing this model in a spreadsheet environment, the annual simulation at an hourly time step took less than one second to complete. The simplicity does come at a price, however, and that price is accuracy. The assumption that the ground temperature near the pipe is not affected by the presence of the pipe goes against what most ground source heat pump developers count on, and that is heat storage in the ground. This model will be tested against more rigorous models to determine how much accuracy is lost in the assumptions and simplifications made. 4.2 Model Comparisons: VBA Earth Tube vs. Fluent The earth tube model, as implemented in VBA, requires the inputs found in table 4.1. To reproduce the model in Fluent, modifications to the Fluent model discussed in chapter 3 are required. Geometrical constraints such as pipe diameter, depth, length, and thickness are all very similar. In the VBA implementation, there is no finite difference solution, so the domain is simply the pipe along with a representative soil temperature with which the pipe interacts. For the Fluent model, however, the geometric constraints must include distances to apply boundary conditions in the 3D domain. These distances represent the width and height of the soil that is modeled. The soil is modeled underneath the pipe, and around the side of the pipe. The VBA earth tube model does not take into account any heat transfer from these regions. In an effort to produce the most direct comparison of the two models, the most 33 appropriate way to implement these boundary conditions in the Fluent model is to use adiabatic boundary conditions on these faces. The soil and fluid properties may be directly placed within both models, which helps ensure agreement between the two. The fluid to pipe convection coefficient will cause some difference between the models. For the earth tube model, a heat transfer correlation is used which is based on Reynolds number, and therefore based on bulk flow rate through the pipe. In Fluent, the local velocity and boundary layer will be simulated, and so a uniform convection coefficient will not be modeled. The last item that changes is the soil temperature prediction. The earth tube model uses annual average surface temperature data to produce a single correlation as a function of depth and time. This then produces a temperature in the soil near the pipe at each time step of the simulation. In Fluent, this would over specify the problem, since the soil is being modeled as a solid volume, not just the near pipe vicinity. The inputs to the earth tube model are the annual average surface temperature and the annual surface temperature amplitude. These two inputs can also be used in Fluent to develop a userdefined transient boundary condition. Some details on userdefined functions are presented in section 3.1. The surface temperature can be approximated as a sine wave throughout the year with the same given average and amplitude. In summary, the models can be made very similar, but with three significant differences: 1. The soil is modeled from the surface to a distance below the pipe and in the region to the side of the pipe in Fluent. The ’deep ground’ and ’farfield’ boundaries will initially be modeled as adiabatic to account for the fact that the VBA earth tube model ignores heat transfer interaction with these regions. 2. The air flow through the tube is modeled as a finite element mesh in order to predict heat transfer behavior from the air in Fluent. In the earth tube model, a correlation is used based on bulk flow through the tube. 34 3. The soil surface boundary condition will be a userdefined sine wave approximation of a yearly cycle with a given average and amplitude in Fluent. In the earth tube model, the same average and amplitude are used in a correlation to predict the temperature of the soil in the nearpipe vicinity, which then becomes the boundary condition for the model. In other words, the earth tube calculations only use the surface temperature to precalculate the soil temperature at the pipe depth for any given day of the year, and the surface is not directly modeled. Note that the VBA implementation was not directly verified against the EnergyPlus implementation. A paper by Lee & Strand (2006) covers the development of the model and implementation of the model in EnergyPlus. The work in this paper was used as the guideline for the current VBA model implementation. One major difference that may occur between the VBA and EnergyPlus implementation is the time step management; the VBA implementation uses a constant one hour time step, while EnergyPlus varies the time step as needed throughout the simulation. The model is updated with the current ground temperature and the earth tube is simulated. Based on a verification of the governing model equations, the model presented in the paper and the VBA model should produce the same results. 4.3 Model Comparison Results: VBA Earth Tube and New Fluent Earth Tube With the proper mesh developed for the Fluent model, the earth tube model as implemented in VBA and the Fluent model can be directly compared. The detailed differences between the models can be found in Section 4.2. The inputs used which are common to both models are given in Table 4.5. The inputs which are not shared between the models are summarized below. For each input, the method used by each model is given. 35 Table 4.5: Common Input Data For Earth Tube Models Input Description Value Units Pipe Data Burial Depth 2.5 [m] Length 25.0 [m] Radius 0.25 [m] Thickness 0.01 [m] Thermal Conductivity 200 W m−K Soil Data Thermal Conductivity 1.5 W m−K Fluid Data Working Fluid Air [−] Flow rate 0.785 m3 s Simulation Data Surf Temp Avg 15 [C] Surf Temp Ampl 5.6 [C] • Ground Surface Temperature – VBA Earth Tube: Use average surface temperature and amplitude data along with equation (4.1) to generate soil temperature near pipe. – Fluent Earth Tube: Use average surface temperature and amplitude data as a transient sine wave: Tsurf (t) = Tsurf,avg + Tsurf,amplsin (A ∗ t). The A coefficient is modified to fit the time unit conventions in Fluent. This equation is similar to the full Kusuda & Achenbach (1965) correlation, but does not require the exponential because the depth is zero at the ground surface. 36 • Deep Ground Heat Transfer – VBA Earth Tube: This model does not take into account any heat transfer interaction to deep ground or heat transfer to the side of the pipe. – Fluent Earth Tube: Soil surrounding pipe is modeled, but all surfaces are set as adiabatic to mimic current earth tube behavior. • Fluid Properties – VBA Earth Tube: This model uses constant values of air properties such as density and specific heat. – Fluent Earth Tube: Fluent contains libraries of material properties, some of which are constant, and some of which are calculated by equations of state (such as ideal gas). – Note: For this model comparison, the air will not see drastic temperature or property changes, so this change is not a significant source of error. The main output of the VBA earth tube model is the outlet temperature from the tube. The earth tube model inherently assumes steady state behavior. During an annual simulation, the ground temperature and inlet temperature are updated, and a new steady state solution is found. For both models, a time step of one hour is used. The simulation is then operated for 8760 time steps to include a full annual simulation. In the earth tube model, hourly weather data is read directly. For the Fluent model, the entering air temperature is a boundary file, which consists of ’ifthen’ structures to select the outdoor dry bulb temperature for any given hour of the year. The weather data used is from Atlanta, GA. Although Fluent can be used to output any data desired, the variable of most interest is the exiting fluid temperature. This is what will be compared to the earth tube model. The goal of this comparison is to determine how well the earth tube 37 Figure 4.1: Earth Tube Models: Exit Air Temperatures can predict ground heat transfer, compared to the more robust and thorough Fluent model. Other comparisons are performed as needed. The exiting air temperatures are plotted for the entire simulation year in Figure 4.1. The numerous simplifications inherent in the earth tube model result in some disagreement with the Fluent model. The deviations between the models are plotted in figure 4.2. The earth tube model is a steady state model at each time step, while Fluent solves the transient equations of momentum and energy for each time step. At the beginning of the simulation, Fluent initializes the entire domain to a steady condition based on the values at the initial time step. For the entire year, the trends of both models agree, but the actual values do not. For most of the year, the models deviate by about 5◦C. In the summer, the Fluent model shows better precooling capabilities, while in the winter, it shows better preheating capabilities. The earth tube model appears to be in error, but in a conservative manner. This seems to be due to the better ground heat storage modeling in Fluent. Another way to characterize the error between the two models is to show the different heat transfer rates predicted. This is similar to presenting the pipe delta 38 Figure 4.2: Earth Tube Models: Exit Air Temperature Deviations temperature (inlet to output), but gives a more applicable aspect to the results. Most designers will care more about the heat transfer rate than predicted temperatures. Figure 4.3 shows the heat transfer rates predicted by the two models. The heat transfer rates between the two models show the same basic trends, but differ drastically in actual values. Through the simulation, the heat transfer rate predicted by Fluent is up to four times as high as the earth tube model prediction. The model predicts higher heat transfer rates in both heating and cooling conditions. This shows that the earth tube model is, in fact, a conservative model. When a designer is using the earth tube model in EnergyPlus within a simulation, it should be noted that the earth tube may actually perform better than the model predicts. Although the earth tube model does deviate from the Fluent model, the extreme increase in simplicity for model implementation and computation speed still make it a possible candidate to become a buried water pipe model. A Fluent model will also be developed to verify the model accuracy when the earth tube is converted to a buried pipe model. This model development is described in section 4.4. 39 Figure 4.3: Earth Tube and Fluent Heat Transfer Rates 4.4 Earth Tube Conversion to Pipe Model This section describes the changes made to the earth tube model in order to simulate a buried pipe with water flow. For simplicity in this study, many simulation parameters will remain constant. The surface temperature will still be modeled by a sinusoidal wave in Fluent, using the same average and amplitude as the earth tube model. The pipe depth will remain constant, although the pipe diameter will change. The new pipe will be a 0.1 [m] nominal diameter pipe. The other boundary conditions will remain constant except for the inlet conditions. For this, a new set of boundary values based on ground source heat pump water temperatures will be used. This data was recorded at the Oklahoma State University Heat Pump Laboratory. The data is the entering water temperature to a set of vertical boreholes, so the data cannot be used as validation for a horizontal heat exchanger loop. However, the data does provide a semirealistic input to the model. There is one key issue with simulating a short section of pipe and not an entire loop. The error in exiting temperature will have a tendency to propogate through the entire loop, and the entering water tem 40 perature will be a function of the exiting water temperature. By using this entering temperature as the boundary condition, the small amount of error in the loop is corrected each new hour. If the full loop were run along with a heat pump model, this error may accumulate. This simulation is out of the context of the Fluent reference model, and requires long term validation data to account for these long term effects. However, as the current model is only for a buried pipe, using this entering water temperature data will ensure that the simulation is performed in the most realistic manner possible without drastic changes to the model itself. The values are summarized in appendix C, and the same values will be used in the Fluent model and the earth tube model. The water flow rate will also be different from the air flow rate. The water velocity will be decreased from 1.0 m/s to 0.15 m/s. For a pipe of 0.1 m diameter, this is on the order of 0.0012 m3/s, or 19 GPM. For typical standard conditions, this results in a Reynolds number of approximately 12,000. For the water flow, Fluent will automatically solve the momentum and boundary layer equations using water data instead of air data once this is selected as the working fluid. For the earth tube model, a new convection correlation will be required. With this Reynolds number known, the convection correlation can be specified in equation (4.8), with the the correlation applicable over the range specified in Table 4.6, as found in Incropera & DeWitt (2002). NuD = 0.027Re4/5 D Pr1/3 μ μs 0.14 (4.8) For the case of this simulation, the minimal property variations will not have a drastic effect on the results because the temperatures throughout the entire simulation vary only within a small range. This assumption should allow the last term in equation (4.8), μ μs 0.14 , to be neglected, thus simplifying the simulation. With this convection correlation in place, the simulation can be performed. For an annual simulation with the given data, the resulting temperatures are shown in figure 41 Table 4.6: Convection Correlation Information for Circular Pipe Flow Conditions: Turbulent Flow Fully Developed 0.7 ≤Pr≤ 16, 700 ReD ≥ 10, 000 L/D≥10 Variable Description ReD Reynolds number based on diameter Pr Fluid Prandtl number μ Viscosity at fluid bulk temperature μs Viscosity at film temperature 4.4. There are three curves shown: Pipe Inlet Temp, Pipe Outlet Temp, and Temp Deviation. The inlet and outlet are basic flow temperatures, although on this plot, the differentials are actually hard to see. With this simulation, the temperature difference from inlet to outlet is small, on the order of 1◦C. This is why the differentials are also plotted. The temperature difference is plotted on a different scale to make it easier to visualize. As can be seen on the plot, the temperatures fluctuate widely. This is simulating the heat pump demand during the week, and then off on the weekends. The model predicts about 1◦C of precooling in the summer, and 0.5◦C of preheating in the winter. The Fluent model was also modified, as discussed previously. The main change in the mesh involves a smaller pipe diameter. This smaller pipe requires a tighter mesh as well, but this is a simple task to perform once a baseline journal file has been created for Gambit. The new pipe mesh was placed in Fluent, and the same simulation that was performed for the earth tube was repeated, with two changes: Pipe inlet temperature and working fluid. The simulation was initialized to data 42 Figure 4.4: Earth Tube Pipe Model Results at the first time step until steady state, and then the transient annual analysis was performed. The exit temperatures for both the Fluent pipe model and the earth tube pipe model are shown in figure 4.5. Figure 4.5 shows that the Fluent model predicts a much higher precooling rate during the summer, and a significant but lower preheating rate during the winter. To better visualize the differentials, the absolute value of the deviations is plotted in figure 4.6. This difference means that, as with the earth tube model comparison, the earth tube is a conservative simulation. If used as a design tool, it appears that the actual performance will be higher. This is most likely due again to the Fluent model capability to better model heat storage in the ground. For the earth tube, the air has an effect on ground temperatures near the pipe, but this effect is small relative to the effect from a water pipe. For this water pipe, the heat storage in the ground is more prominent. Improper modeling of the ground capacity along with previously discussed model assumptions led to poor results from the modified earth tube model. 43 Figure 4.5: Fluent and Earth Tube Pipe Model Exit Temperatures Figure 4.6: Fluent and Earth Tube Pipe Model Delta Exit Temperatures 44 4.5 Earth Tube Model Guidance The results of the pipe study show that the earth tube model does not give good results when using water as the working fluid. One reason for this is the added thermal storage in the system. In a typical earth tube, the air flows through quickly, and the effects of the air are not as prominent in the ground as a flow of water. The model creates a steady state simulation at each time step, assuming a constant temperature in the ground near the pipe. With water, the ground absorbs and rejects significant heat which causes the variation in temperature along the length of pipe to be more intensified. An improvement to this model could be found in the calculation of the ground temperature in the near pipe region. Currently, surface temperature data is used in a correlation and this results in a temperature at the pipe depth, which is used as a boundary condition to the model. This temperature, however, is too close to be considered a farfield boundary. A serial model is not robust enough, and there should be some feedback from the pipe when generating a ground temperature. This is a difficult improvement when the model remains steady state. The earth tube model is built upon deep underlying assumptions. Although the model has limitations, some general guidance is provided for situations where the earth tube model should provide higher accuracy. Due to the constant temperature assumption along the length of the pipe, the earth tube model accuracy should be improved in designs where the earth tube will not encounter major temperature variations through the ground along the length of the tube. This is similar to the statement that the earth tube model will be more accurate for cases where ground temperatures are dominated by groundsurface interaction rather than groundpipe interaction. • Minimal earth tube pipe length: The shorter the tube, the less chance there is for a temperature variation within the tube, and therefore within the ground. 45 When a shorter tube is in place, however, the entrance effects are much more definitive. The heat transfer coefficient will be higher in this entrance region due to the developing boundary layer, but overall, the heat transfer will be reduced. • Increased flow rate: When the flow rate is reduced, the temperature difference in the fluid will be greater. This axial variation poses an opportunity for inaccuracy in the model. Therefore, it may be inferred that the model will increase accuracy with higher flow rate. • Soil properties: With the earth tube model assuming steady state behavior at every time step, higher accuracy should be found in simulations that better approximate steady behavior. Higher values for soil conductivity and diffusivity lead to faster heat diffusion through the soil. This will in turn reduce the possibility for temperature variation in the soil. This will become a case that performs similar to a steady situation. With regard to soil properties, the higher diffusivity used should result in the higher accuracy model. In addition to the model input guidance given above, a parametric study was performed with the earth tube model within the EnergyPlus simulation program to determine which of the model inputs have a significant impact to the simulation output. For the inputs with higher sensitivity, more time should be spent increasing the quality. If an input shows negligible sensitivity, time will be wasted trying to improve its quality. A simple earth tube model was set up for the annual simulation. The simulation contains a full building simulation, and the earth tube is applied to a single zone. The representative output for the model is selected as the furnace heat transfer rate. One at a time, the variables shown in table 4.7 were adjusted, and the change in coil energy requirements were recorded. The sensitivity was nondimensionalized according to equation (4.9) with inputs X and outputs Y. The subscripts refer to 46 either 0 for the base case simulation or i for the simulation of case i. Sensitivity = Yi−Y0 Y0 Xi−X0 X0 (4.9) With equation 4.9, the study revealed the data in table 4.7. Table 4.7: Earth Tube Model Sensitivity Study Results Case Variable Baseline New Inp Input Chg Output Chg Sensitivity 0 Baseline      1 Flow (m3/s) 1 1.25 0.25 0.0007 0.0028 2 Radius (m) 0.3048 0.4048 0.3281 0.0003 0.0008 3 Length (m) 50 75 0.5 0.0063 0.0127 4 Depth (m) 3 4.5 0.5 0.0002 0.0003 In addition to testing the sensitivity of the model, influence coefficients can provide another level of meaning to how the model responds to different inputs. Application of influence coefficients to building simulation are discussed in Spitler et al. (1989). These coefficients are similar to the sensitivity discussed previously. However, the estimated percent error is developed which includes the estimated error (uncertainty) in the input value. Equation (4.10) is used to define the type 2 dimensional influence coefficients, and equation (4.11) is used to define the percent error. Influence Coefficient = Y0−Y1 Y0 X0 − X1 (4.10) Percent Error = Influence Coefficient · Estimated Error in Input · 100 (4.11) The influence coefficients defined here are dimensional as they contain the inverse units of the input parameter. When the percent error is calculated, the estimated 47 error in input carries the input units also, which results in a dimensionless percent error, as expected. The four case results are presented in Table 4.8 in terms of influence coefficients and percent error. The uncertainty of the inputs was set to 5%. Table 4.8: Earth Tube Model Influence Coefficients and Percent Error Case Variable Influence Coefficient Percent Error 1 Flow (m3/s) 0.00276 1.378% 2 Radius (m) 0.00251 1.253% 3 Length (m) 0.00025 0.127% 4 Depth (m) 0.00012 0.058% This influence coefficient study shows that the percent error in results is highest in the flow rate and pipe radius. The sensitivity of the model was shown above to be higher in the flow rate and pipe length input parameters. With this information, it is assumed that the highest possibility for model improvement is in the flow rate estimation. The simulation indicates that the earth tube model is sensitive to the flow rate, and a small uncertainty in the input parameter produces a percent error of about 1.4%. 4.6 Earth Tube Model Conclusions The results of this chapter indicate that the earth tube model, as it operates currently, is not a suitable candidate for the buried pipe model. Other models will be compared throughout this study, most of which include a more substantial set of input data and require significantly more computation time. If none of these seem suitable, the earth tube model may become the basis of a new model that better handles heat storage in the soil and pipe. 48 CHAPTER 5 Mei Model 5.1 Model Description Mei (1986a) developed a model to predict heat transfer around a buried horizontal pipe. This model uses an explicit finite difference scheme, and develops update equations, which are used in a pseudo3D environment along the pipe. As one would expect, this model uses symmetry to split the soil region in half, centered along the pipe. This model is considered simple due to a single cylindrical coordinate system. This gives a less complex algorithm within the code compared to other models which may mix coordinate systems. The original code, however, seems more complex, as the model was programmed using older FORTRAN conventions. This model was translated to VBA and tested in the Excel environment. There are actually numerous advantages when using this model. One, already mentioned, is the single coordinate system. This simplifies the heat transfer equations in the domain. Another advantage this model has is that it allows for a backfill material to be used around the buried pipe. This ability is already in the code, so the only items a user must input are the material properties of the backfill material, and how thick the material is, and the program will automatically simulate the region accordingly. In this same manner, if no backfill material is used, this region can simply be set to the outer soil properties, and the program will run as normal, with soil filling the backfill region. There are also some disadvantages of this model. The single radial coordinate system is centered on the centerline of the pipe, and extends the radius to the ground 49 surface (although the domain can be fully underground). The model layout can be seen in figure 5.1. This layout indicates that there will be a single node at the ground surface at each cross section. This does not seem to be as effective as spreading out nodes along the ground surface. However, testing may indicate that a single node is suitable for representation. Another disadvantage may be the distribution of grid points. This model does allow the grid points to be nonuniform in the radial direction in the soil, with fixed spacing in each of the backfill and pipe wall regions. However the variation in grid spacing can lead to stability problems. The radial coordinate system is applied to nodes in different ways throughout the domain. The water at any given pipe cross section is represented as a single node. The pipe wall is represented by radial nodes; variations in the angular (θ) direction are ignored. This correlates to the assumption that the heat transfer up to the outer pipe wall is axisymmetric. Therefore, at a given cross section and radius, the pipe temperature is assumed constant for all angles. Moving outward radially, the pipe outer wall interfaces with the backfill region. The code has an interesting effect through this backfill region. There are nodes set up as if a full 2D (radial and angular) grid was set up, but the heat transfer only passes in the radial direction. This is analogous to breaking the domain up into a 3D grid, but ignoring the heat transfer along the axial direction. This allows a firststage analysis of heat transfer varying with angle and therefore depth. Moving outward the next section is the main soil section. In this region, the domain is fully 2D, both angular and radial heat transfer effects are modeled. To better illustrate this grid, figure 5.1 describes the various domains and figure 5.2 shows the layout of the finite difference grid. As mentioned previously, there are a number of assumptions that are used by the model, in order to develop a tractable solution. Some of these may turn out to be invalid in certain situations. The main assumptions used, as presented by Mei, are summarized as follows: 50 Figure 5.1: Mei Model Domain Description Figure 5.2: Mei Model Finite Difference Grid 51 • Homogeneous soil, which allows constant properties within the soil region. However, with the addition of the backfill region, a secondary homogeneous material (either a second soil type or pipe insulation) can be supplied for the analysis. • Fluid temperature and velocity are constant at any cross section. This allows a single node to represent the fluid temperature and velocity within the pipe at any axial point. • Pipe depth is great enough that the domain boundary can be considered a farfield boundary. This assumption allows the model to ignore surface heat transfer effects directly. • There is only one pipe in the ground in the given domain. • Heat transfer within the pipe up to the pipe wall is axisymmetric. This allows the heat transfer within this region to be 1D radial. The model operates on a dual time step. The pipe wall and water nodes run at a shorter time step in an attempt to ensure convergence, while the backfill and soil nodes run at a much longer time step. This allows the accuracy and stability of a short time step only in the regions where it is required. In the backfill and soil regions, where the temperature gradients are smaller, a longer time step is more suitable. As one would expect, since the model is broken down into four different regions (water, pipe wall, backfill, and soil), there are four different temperature arrays which hold the data. The water temperature array is defined in the program as U (Length, Time Step). The pipe wall temperature array is defined as T (Radius, Length, Time). In the backfill region, the heat transfer is modeled as pseudo 2D; the temperature can vary angularly, but the heat transfer along the angular direction is ignored. The array will therefore have four variables, and is defined as F (Radius, Length, Angle, Time). The final temperature array is the soil temperature 52 Figure 5.3: Mei Model Update Equation Locations array. It is defined as S (Radius, Length, Angle, Time). These arrays are updated sequentially, according to specific update equations for the simulation domain and model conditions. The flow of the model in its entirety is similar to other models. The model reads parameters supplied by the user, initializes the domain, starts the transient iterations, updates boundary temperatures, updates the water and pipe wall temperatures on the smaller time step and then updates the soil and backfill regions on the larger time step, writes output at each larger time step, and repeats until the simulation is complete. More detail is given later in this section as required. The update equations differ based on the region and conditions. For instance, the pipe wall to water interface equations change depending on whether the flow is on or off. The different locations to which an update equation might be applied are shown in figure 5.3. Table 5.1 references the locations shown in figure 5.3, describes the conditions, and references the equation number of the update equation discussed in the following paragraphs. In the equations, the independent variables have been simplified. The time step is not listed as a subscript as are the other variables. Instead, the convention has been 53 Table 5.1: Mei Model Update Equation Organization Location Number Description (Figure 5.3) (in text) (Conditions) A 5.1 Water Node, Flow ON, J=1 A 5.4 Water Node, Flow OFF, All Other J B 5.3 Inner Pipe Wall, Flow ON B 5.5 Inner Pipe Wall, Flow OFF C 5.6 Internal Pipe Wall Node, All Conditions D 5.7 Outer Pipe Wall, All Conditions D 5.8 Inner Backfill Node, All Conditions E 5.9 Internal Backfill, All Conditions F 5.10 Outer Backfill, θ = 0 F 5.11 Outer Backfill, 0 < θ < π F 5.12 Outer Backfill, θ = π G 5.13 Soil, θ = 0 G 5.14 Soil, 0 < θ < π G 5.15 Soil, θ = π H 4.1 Farfield, All Conditions 54 Table 5.2: Mei Model Nomenclature Parameters and Subscripts Region Extents Pipe Length Inlet to Outlet J 1 to NZ Angular Location Top to Bottom of Domain M 1 to NRAD Pipe Wall Radial Location Pipe Centerline to Outer Pipe Wall I 1 to NPIPE Backfill Radial Location Outer Pipe Wall to Outer Backfill I 1 to NFZ Soil Radial Location Outer Backfill to Outer Domain Boundary I 1 to NSOIL Example Variable Description UJ Water Temp at Length J TI,J Pipe Wall Temp at Radius I, Length J FI,J,M Backfill Temp at Rad I, Length J, Angle M SI,J,M Soil Temp at Radius I, Length J, Angle M applied which uses the prime ’ notation to denote the updated time step. Aside from the time step, the independent variables are given in table 5.2. As mentioned, the variables in this section represent both the updated temperatures and the temperature from the previous time step. The updated temperatures are denoted with a prime, such as F I,J,M. Note that with this dual time step program, the main outer time step, which runs the soil and backfill regions, is denoted Δt1. The inner time step, which runs the coil wall and water nodes, is denoted Δt2. These variables will appear throughout these equations. In this domain, there are a number of radii which are very specific. The general 55 radius term is denoted rI , and can be used within any of the regions. There are also four specific radii that represent specific interfaces in the domain. These, along with the last variables of interest in these equations, are denoted and defined in the following: • r1: Inner Pipe Wall Radius • r2: Outer Pipe Wall Radius / Inner Backfill Radius • r3: Outer Backfill Radius / Inner Soil Radius • r4: Outer Soil Radius / Farfield Boundary Radius • h: Water to Pipe Convection Coefficient • V : Flow Velocity • Δr: Radial spacing: This variable has a constant value in the pipe wall and backfill regions, but is nonuniform in the soil region. Therefore, the radial spacing is dependent on the current radius. The first set of equations presented here are applied to the water and inner pipe wall nodes when the flow is on. U 1 = 1 − V dT dz − 2hΔt2 ρfCpf r1 U1 + V dT dz E.W.T. + 2hΔt2 ρfCpf r1 T1,1 (5.1) U J = 1 − VdT dz − 2hΔt2 ρfCpf r1 UJ + V dT dz UJ1 + 2hΔt2 ρfCpf r1 T1,J (5.2) T 1,J = 1 − αpipeΔt2 Δrpipe 2 Δrpipe + 1 r1 − 2αpipeΔt2h Δrpipekpipe T1,J + αpipeΔt2 Δrpipe 2 Δrpipe + 1 r1 T2,J + 2αpipeΔt2h Δrpipekpipe UJ (5.3) 56 Equations (5.4) and (5.5) are used when the flow in the pipe is off. U J = UJ + 2kpipe (T2,J − T1,J)Δt2 ρfCpf r2 1log 1 + Δrpipe r1 (5.4) T 1,J = U J (5.5) Once the simulation has left the inner pipe wall surface, the equations are no longer dependent on the flow being on or off. Therefore, all of the remaining equations in this discussion are applied for both flow on and flow off conditions. Equations (5.6) and (5.7) are for the rest of the pipe wall nodes, both interior nodes and outer pipe wall surface node, respectively. In equation (5.7), a term TAVG has been introduced. Within the pipe wall, there is a single set of pure radial nodes. Outside the pipe wall, there are multiple nodes in the angular direction, even though angular interaction is neglected. Therefore at the outer pipe wall interface, the grid turns from a single 1D radial grid to a number of radial grid points. The single interior node must interact with the multiple grid points outside the pipe wall. Therefore, it is assumed that the one outer pipe wall node interacts with the average of all innermost backfill nodes. This average of the backfill nodes which are directly adjacent to the pipe wall is called TAVG for brevity. T I,J = 1 − αpipeΔt2 Δrpipe 1 Δrpipe + 1 rI − αpipeΔt2 Δr2 pipe TI,J + αpipeΔt2 Δrpipe 1 Δrpipe + 1 rI TI+1,J + αpipeΔt2 Δr2 pipe TI−1,J (5.6) T Npipe,J = 1 − 2 kfill kpipe αpipeΔt2 ΔrpipeΔrfill − αpipeΔt2 Δrpipe 2 Δrpipe − 1 r2 × TNpipe,J + 2 kfill kpipe αpipeΔt2 ΔrpipeΔrfill TAVG + αpipeΔt2 Δrpipe 2 Δrpipe − 1 r2 TNpipe1,J (5.7) 57 When the simulation has reached this point, both the water and the pipe wall nodes have been simulated. This occurs at a time step of Δt2. The simulation then moves on to the backfill region and simulates at a time step of Δt1. There are three subregions in the backfill region, the inner node, interior nodes, and the outer node. Each of these are defined individually in equations (5.8) through (5.12). F 1,J,M = T J,NPipe (5.8) F I,J,M = 1 − αfillΔt1 Δrfill 1 Δrfill + 1 rI − αfillΔt1 Δrfill 2 FI,J,M + αfillΔt1 Δrfill 1 Δrfill + 1 rI FI+1,J,M + αfillΔt1 Δrfill 2 FI−1,J,M (5.9) Once the simulation has reached the outermost backfill node, the grid has three possible equations based on the node location. Even though the update equations in the interior of the backfill region are only derived in the radial direction, the soil region is derived in both the radial and the angular directions. The grid points along the interface of the backfill and soil regions are derived to match the soil region. Because of the angular dependence, the grid points at the angular domain extents will encounter the adiabatic (symmetric) boundary. This leads to three different update equations, one for M=1 (θ = 0 radians at the top of the domain), one for M=NRAD (θ = π radians at the bottom of the domain), and one for the interior points. These locations are described in figure 5.4. The equations follow. 58 Figure 5.4: Mei Model Description: Angular Grid Locations F NFZ,J,1 = 2ksoilαfillΔt1 kfillΔrsoil,1Δrfill S2,J,1 + 2αfillΔt1 Δrfill 1 Δrfill − 1 2r3 FNFZ−1,J,1 + αfillΔt1 (r3Δθ)2 2FNFZ,J,M+1 + FNFZ,J,1 × 1 − 2ksoilαfillΔt1 kfillΔrsoil,1Δrfill − 2αfillΔt1 Δrfill 1 Δrfill − 1 2r3 − 2 αfillΔt1 (r3Δθ)2 (5.10) 59 F NFZ,J,M = 2ksoilαfillΔt1 kfillΔrsoil,1Δrfill S2,J,M + 2αfillΔt1 Δrfill 1 Δrfill − 1 2r3 FNFZ−1,J,M + αfillΔt1 (r3Δθ)2 (FNFZ,J,M+1 + FNFZ,J,M−1) + FNFZ,J,M × 1 − 2ksoilαfillΔt1 kfillΔrsoil,1Δrfill − 2αfillΔt1 Δrfill 1 Δrfill − 1 2r3 − 2 αfillΔt1 (r3Δθ)2 (5.11) F NFZ,J,NRAD = 2ksoilαfillΔt1 kfillΔrsoil,1Δrfill S2,J,M + 2αfillΔt1 Δrfill 1 Δrfill − 1 2r3 FNFZ−1,J,M + αfillΔt1 (r3Δθ)2 2FNFZ,J,M−1 + FNFZ,J,NRAD × 1 − 2ksoilαfillΔt1 kfillΔrsoil,1Δrfill − 2αfillΔt1 Δrfill 1 Δrfill − 1 2r3 − 2 αfillΔt1 (r3Δθ)2 (5.12) With these equations applied at M=1 through M=NRAD, the backfill region outer interface is complete, and the simulation can move to the soil region. The soil domain is meshed into a full 2D polar grid. This mesh is similar to the interface of the backfill and the soil in that there will be three different update equations based on the angular location. These three equations account for the adiabatic surfaces at the top and bottom of the left side of the domain, as well as the internal points. This mesh structure can be visualized as before in figure 5.4. Equations (5.13), (5.14), and (5.15) are for M=1 points, Internal points, and M=NRAD points, respectively. 60 S I,J,1 = 1 − 2αsoilΔt1 Δrsoil,IΔrsoil,I1 − αsoilΔt1 Δrsoil,I1 2 Δrsoil,I1 − 1 rI − αsoilΔt1 (rIΔθ)2 × SI,J,1 + 2αsoilΔt1 Δrsoil,IΔrsoil,I1 SI+1,J,1 + αsoilΔt1 Δrsoil,I1 2 Δrsoil,I1 − 1 rI SI−1,J,1 + αsoilΔt1 (rIΔθ)2 2SI,J,2 (5.13) S I,J,M = 1 − 2αsoilΔt1 Δrsoil,IΔrsoil,I1 − αsoilΔt1 Δrsoil,I1 2 Δrsoil,I1 − 1 rI − αsoilΔt1 (rIΔθ)2 × SI,J,M + 2αsoilΔt1 Δrsoil,IΔrsoil,I1 SI+1,J,M + αsoilΔt1 Δrsoil,I1 2 Δrsoil,I1 − 1 rI SI−1,J,M + αsoilΔt1 (rIΔθ)2 (SI,J,M+1 + SI,J,M−1) (5.14) S I,J,NRAD = 1 − 2αsoilΔt1 Δrsoil,IΔrsoil,I1 − αsoilΔt1 Δrsoil,I1 2 Δrsoil,I1 − 1 rI − αsoilΔt1 (rIΔθ)2 × SI,J,NRAD + 2αsoilΔt1 Δrsoil,IΔrsoil,I1 SI+1,J,NRAD + αsoilΔt1 Δrsoil,I1 2 Δrsoil,I1 − 1 rI SI−1,J,NRAD + αsoilΔt1 (rIΔθ)2 2SI,J,NRAD−1 (5.15) Thus all interior points of the calculation domain have been acknowledged. The only points which do not have an update equation given here are along the outer farfield boundary. The depth of each boundary point in this radial domain is calculated 61 based on input geometry. With the current node depth and simulation time known, the Kusuda and Achenbach correlation which has been discussed previously in this work is used to update the node at each outer time step. The Kusuda correlation used here is very similar to equation (4.1). More detail is provided in section 5.3.4. 5.2 Modifications to Mei Model After the initial implementation of the Mei model, some modifications were required to fit better with the current implementation language and some features were removed or modified. The structure of the original program was very difficult to understand because the programming did not include some common structures such as the ’Else If’ command. Without a structure such as this, the code must wind in and out with hundreds of ’GOTO’ statements, rendering it quite difficult to follow. Thus the program was modified and updated with newer command structures whenever possible. As for actual model features, there were two items which were modified: specification of the backfill region, and temperature arrays to store previous values of backfill and soil nodal temperatures. In the original program, the use of a backfill region was optional. This allowed a user to place the pipe in direct contact with the soil. However, this option introduces many additional equations. In order to simplify the program, the option was removed, and the backfill region is now in place all the time. This does not actually limit the program in any way. To simulate the pipe without a backfill region, one can simply specify the backfill properties to be the same as the soil properties. The second change involved the storage of values from previous time steps. Although it may have been a simple omission, the code in Mei’s report never included storing the backfill and soil temperatures at the end of each time step. This was revealed during testing as discussed in section 5.3.3. The addition of storage resulted in a model that performed as expected. 62 5.3 Model Testing After initially debugging the model, it appeared to be operating correctly. However, tests were performed to help reveal any underlying model problems. The first test was a simple test where the entering water temperature and all boundary conditions were set to the same temperature. If the model were operating correctly, there would be essentially no response from the model, just constant temperatures everywhere. For the second test, the inputs to the Mei model were adjusted to mimic a constant pipe wall temperature condition. This was set by carefully adjusting the conductivities and other parameters. In the third test, the model was run for different initial conditions to determine convergence characteristics and to ensure that the model was properly updating each temperature node. In the fourth test, the flow rate was increased to eliminate axial heat transfer effects and to assess the 2D effects calculated by the models. 5.3.1 Fully Constant Model Testing This test was designed to ensure that the model algorithm and update equations are defined properly. If the update equations are not implemented properly, heat transfer may be simulated through nodes when it should not, simply because of an error in a coefficient. During this test, the boundary temperatures and inlet water temperature were set to a constant value of 50◦F. The model responded properly to these boundary conditions. This test served to establish some confidence in the model implementation. 5.3.2 Constant Wall Temperature Testing This test approximates an analytic solution. The analytic solution is simply a constant pipe wall temperature solution for exiting fluid temperature as a function of entering water temperature. A more detailed discussion of the analytic solution is 63 given in section 6.4.2. In this test, the Mei model was modified to approximate a constant pipe wall. The conductivities in the backfill region, soil region, and pipe wall were increased to the order of 999 Btu/hrftF, and the outer boundary temperature was held constant. With the high conductivities, this boundary temperature should transfer quickly to the inner pipe wall. This test was performed with a boundary temperature of 50◦F and an entering water temperature of 0◦F. The test was repeated for several different pipe lengths to see the model behavior. At each run, the pipe inner surface temperature along the length of the pipe was recorded. With this numerical scheme, it is nearly impossible to get the pipe temperatures to exactly 50◦ while achieving model convergence. Therefore, an average surface temperature for each run was found, and this was used as the surface temperature in the analytic solution. The exiting water temperature for the numerical case was then compared to the analytic solution. The main analytic equation governing the exiting fluid temperature is given here as equation (5.16), but as mentioned, the full development is in section 6.4.2. As is also mentioned in that section, the analytic solution was applied to small segments of the pipe in order to improve the comparison with the Mei model. T2 = k0.023 ρDV μ 0.8 Pr0.4πL Ts − T1 2 + ˙ mCpT1 k 20.023 ρDV μ 0.8 Pr0.4πL+ ˙ mCp (5.16) The results of this study for the Mei model are given in figure 5.5. At the lengths reported, the models agree to within 5%. This test removes the soil modeling from the simulation, and with this high agreement, increases the confidence in the water modeling of the simulation. The dual time step, which allows the soil to run at a long time step while the pipe wall and water node can run at a short time step, is required for convergence. With a thin pipe, the pipe node spacing becomes very small. When these small nodes are 64 Figure 5.5: Mei Model Test Results: Constant Wall Temperature 65 encountered, a substantially short time step must also be used. During this specific test of constant wall temperature, conductivity values were set above realistic values. This was required in order to quickly translate the boundary temperatures to the pipe wall. The conductivity plays a role in the update equation coefficients, and when these values were increased, the model had serious convergence problems. In order to achieve convergence, the model was run at a 30 minute outer time step, and a .001 minute inner time step. For normal conductivities, initial tests revealed that for an outer time step of one hour, an inner time step of 0.1 minute was required for convergence. The minimum system time step in EnergyPlus is one minute. In order for this model to be implemented in EnergyPlus, an outer structure would be required to simulate this model a number of times during a single EnergyPlus time step. This stability limit does not completely eliminate the possibility of implementing the Mei model in EnergyPlus, but it does introduce another level of complexity to the model and the implementation. 5.3.3 Initial Condition Independence Testing In order to ensure that the model was properly updating the temperatures at each node, a series of initial conditions were placed on the model. The model was then run for a period until the results converged. Figure 5.6 shows the results. In the first test, an initial condition of 0◦F was placed on every node in the model domain. In the second test, the entire model domain was initialized to the Kusuda and Achenbach predicted temperature at the pipe burial depth. In the third test, the model domain was initialized to 100◦F. Figure 5.6 shows that even with the extreme initializations of tests one and three, the model converges, though over two weeks worth of simulation time was required! 66 Figure 5.6: Mei Model Test Results: Initial Condition Independence 67 5.3.4 High Flow Rate Testing (2D Model Simulation) The Mei model does not account for the axial heat transfer effects in the soil or pipe wall. The water in the pipe is the only region which carries heat along the length of the pipe. The Fluent model, however, includes heat transfer in all dimensions. Because of this difference, it is necessary to perform a test which can simulate the ability of the models to simulate in just two dimensions (eliminating the axial dependence). This is performed by substantially increasing the flow rate through the pipe. With this higher flow rate, the fluid will experience almost no temperature change along the entire length of the pipe. If the boundary conditions also do not change along the length of the pipe, two dimensional heat transfer will result. This test marks the first comparison between the Mei model and the Fluent verification model. The Mei model is a cylindrical domain. The Fluent verification model is a rectangular based domain. For many of the boundary conditions, this will not have an effect. In the farfield boundary condition, however, a deviation will occur between the two models. Both models will use the Kusuda and Achenbach correlation to develop the boundary temperatures. In order to create a pseudo steadystate model, the simulation time will be set to zero so that this Kusuda temperature field will not change throughout the iterations. The Fluent model will use the Kusuda and Achenbach directly on the ground surface (no outdoor air interaction), and the farfield (side wall) boundary, and the deep ground boundary. The Mei model will use the Kusuda and Achenbach correlation directly on the cylindrical perimeter of the model. The Mei model does not differentiate between the deep ground and side boundaries, because it is one outer surface boundary in cylindrical coordinates. Figure 5.7 shows the boundaries for both models. While increasing the flow rate, a value will be reached where the axial effects are 68 Figure 5.7: Mei and Fluent Model Boundary Differences negligible. In the Mei model, the fluid enters at 37◦C, with a flow rate of 0.24 m3 s . The flow was increased similarly in the Fluent model. The removal of axial heat transfer from the model should help to reveal two things: whether the models agree when only 2D heat transfer effects are concerned, and if the different boundary type seem to have an effect on overall model agreement. Both of these models are transient in nature, and therefore include the effects of heat storage in the soil and other regions where appropriate. In order to perform this analysis, the models were run until a steady condition was achieved. In the Fluent model, one could turn the model to a steady behavior by switching model parameters, and letting it run until convergence. However, in the interest of comparing these two models in their base configurations, they were both left as transient. To monitor for steady conditions, some variables had to be recorded. For the Mei model, representative temperatures at numerous locations throughout the domain were printed at each time step. In the Fluent model, temperature output is not as convenient, so the overall heat transfer from the pipe wall was recorded, and monitored for steady 69 behavior. Figure 5.8: High Flow Test: Mei and Fluent Cross Section Temperature Profiles Figure 5.8 shows the temperature distributions for both models. As the figure shows, the temperature distributions do not agree well. The temperature distribution along the farfield boundary follows a behavior that would be expected from the Kusuda and Achenbach correlation, but the rest of the domain shows deviations. The Mei model has already shown problems with initial condition independence and stability. The inaccuracy of the 2D calculation is another serious limitation of the Mei model. 5.4 Mei Model Conclusions The Mei model was introduced as a model with intermediate complexity. The single cylindrical coordinate system indicated a model which could produce quality results while balancing the computation time that may be encountered with a more complex model. However, throughout the initial testing, the model did not produce qual 70 ity results and encountered several problems which eliminate it as a candidate for implementation within EnergyPlus: • Stability • Accuracy • Ground Surface Interaction • Initial Condition Independence The first problem, stability, is a result of two things: an explicitly formulated model, and widely varying grid spacing. This varying grid spacing is partially accounted for by the dual time step method used. However, in the initial testing there were several simulations that were unusable because instability caused divergence. This is the main problem of the Mei model. The second problem involves the accuracy of the model. This first set of initial tests, where the domain is idealized in one form or another, show that the Mei model does not agree well with the verification model. The results are expected to worsen as the domain is placed back to its default state where the boundary conditions are both spatially and temporally dependent. The third problem is that the model does not directly interact with the ground surface. The farfield boundary condition of the model is simply a correlation of depth and time. The problem with this is that the model never takes into account the short time scale behavior of the outdoor air. With the other problems encountered during testing, a solution for this problem was not investigated. The fourth problem is that the model showed a very slow independence from initial conditions. The model took almost fourteen days to ’forget’ the initialization of the model. A ’warmup’ period would be required for the model in addition to the base warmup periods that EnergyPlus uses. These problems indicate that theMei model is not currently the ideal model for implementation into EnergyPlus. In chapter 6, the Piechowski model will be evaluated, 71 and depending on these results, chapter 7 will reveal which model is implemented into EnergyPlus. Also, if no model is suitable in its current state, a decision will be made regarding major modifications to one of the models in order to improve the quality, so that it can then be implemented. 72 CHAPTER 6 Piechowski Model 6.1 Model Description Piechowski (1996) developed a quasi3D finite difference buried pipe heat transfer model. Although a set of 3D grid points is generated, the effects of heat transfer from one section of pipe length to the next is only modeled within the fluid. In other words, the soil and pipe material do not model axial heat transfer. Piechowski assumes that the larger heat transfer gradients occur in a radial direction away from the pipe surface, with relatively smaller gradients in the axial direction. Therefore, at each pipe cross section, a 2D grid is developed. The grid is based on a cartesian coordinate system, with boundaries at the ground surface, a deep ground depth, and along the farfield surface on the side. The domain splits the pipe along the middle of the cross section in order to accommodate symmetry. Within this cartesian grid, there is one singular node at each cross section that includes the soil near the pipe, the pipe, and the fluid within the pipe. This layout is shown in figure 6.1. Using a single node alone to represent the pipe would not provide sufficient accuracy, and so a secondary coordinate system was developed, centered upon this representative node. This secondary domain is in radial coordinates. In this radial system, where the grid spacing is more refined than in the outer cartesian system, both heat and mass transfer are calculated, while in the larger grid, only sensible heat transfer is accounted for. The smaller grid layout is shown in figure 6.2. Within the program, there are three distinct temperature arrays. Each of these accounts for a different section of the simulation domain. The first temperature array 73 Figure 6.1: Large Scale View of Piechowski Domain Cutaway Figure 6.2: Small Scale View of Piechowski Domain Cutaway 74 is T (TimeStep,Pipe Length,Soil Depth,Soil Width). This temperature array has a value for the entire Cartesian soil domain. In this region, the grid is a 3D Cartesian coordinate system, and includes a representative value for the rectangular node which contains the pipe. This node will be discussed further in a later section. The second temperature array is held in the variable TR (TimeStep,Pipe Length,Radius). This temperature array represents the values within the radial nearpipe region. In this section, the heat transfer is assumed to be axisymmetric, and the grid is therefore 1D (radial) throughout each cross section. When this region is analyzed in the program, the results are converted into an effective temperature, which can then be used in the Cartesian farfield analysis. In this way, the pipe effects are ’hidden’ from the main routine, which simplifies the entire analysis. The third temperature array is TW (TimeStep,Pipe Length). This array holds all the water temperatures along the pipe length at each time step. The program follows a simple process all controlled by the MAIN subroutine: • Calculate coefficients: pipe convection, moisture coefficients, soil properties • Initialize arrays: water temperature, radial nodes, Cartesian nodes • Loop through all simulation days and perform the following tasks: – Update side boundary soil temperature as a function of depth, time, etc. – Update surface temperature and deep ground surface temperature. – Retrieve cycle ontime for current period in order to provide a stopping point for the simulation. – Retrieve entering pipe water temperature. – Call routine to simulate nearpipe radial system and produce an ’effective’ pipe node temperature for the farfield Cartesian system. 75 – Call routine to develop equation coefficients based on correlations and input data. – Loop through all pipe length cross sections and perform the following tasks: ∗ Loop through all nodes at the current cross section and update all nodal temperatures based on different conditions: Adiabatic on symmetric centerline, constant temperature boundaries on surface, side and deep ground This process is performed with the following required inputs: Material properties, entering water temperatures, pipe oncycle time per day, geometric data, and correlations. Some of the other key features of this simulation from a programming standpoint are the use of a function that can update a boundary temperature for both adiabatic and constant temperature boundaries. Using the same function simplifies the programming, and cleans up the code. For the entering water temperature, the same hourly boundary data was used as in the earth tube model, and is summarized in Appendix C. In Piechowski’s full program, the effects of moisture are modeled in the radial grid system. For this investigation, the effects were removed as described in section 6.2. The heart of the program is the set of algorithms used for updating the temperatures of each node at every time step. The program is designed so that at each tim
Click tabs to swap between content that is broken into logical sections.
Rating  
Title  Development, Verification, and Implementation Of a Horizontal Buried Pipe Ground Heat Transfer Model in Energyplus 
Date  20080501 
Author  Lee, Edwin 
Department  Mechanical Engineering 
Document Type  
Full Text Type  Open Access 
Note  Thesis 
Rights  © Oklahoma Agricultural and Mechanical Board of Regents 
Transcript  DEVELOPMENT, VERIFICATION, AND IMPLEMENTATION OF A HORIZONTAL BURIED PIPE GROUND HEAT TRANSFER MODEL IN ENERGYPLUS By EDWIN LEE Bachelor of Science Oklahoma State University Stillwater, Oklahoma 2006 Submitted to the Faculty of the Graduate College of Oklahoma State University in partial fulfillment of the requirements for the Degree of MASTER OF SCIENCE May, 2008 COPYRIGHT c By EDWIN LEE May, 2008 DEVELOPMENT, VERIFICATION, AND IMPLEMENTATION OF A HORIZONTAL BURIED PIPE GROUND HEAT TRANSFER MODEL IN ENERGYPLUS Thesis Approved: Dr. Daniel Fisher Thesis Advisor Dr. Jeffrey Spitler Dr. David Lilley Dr. A. Gordon Emslie Dean of the Graduate College iii ACKNOWLEDGMENTS There are many people who have helped me along the way. Some have directly helped me with my studies or by providing a bit of humor to ease stressful times. Some have provided inspiration to keep moving in the right direction. It is not conceivable that I could write such a brief summary that would express all the gratitude that I hold for all the ways I have been helped, but I’ll give it the old college try. My utmost praise and worship go to our Lord God Almighty. Without him, there would be nothing, and through him my life has been realized. My mother and father raised me to always be grateful for what I have, and always remember that I have brothers in this world, and guardian angels in the next world watching over me. They have always supported my ventures, even when they seemed sudden and unclear. My brother and sister have always shown me dedication which is unmovable. Throughout tough times, they persevered and created a wonderful family that I am very proud of. They show me understanding when my life seems so busy, and they remind me that there is always time for family. My nieces are seemingly proof that energy is not conserved. They are full of jokes, laughter, smiles, and excitement twentyfour hours a day. Every time I see them, I remember that there is so much in this world to be thankful for, and to never take time for granted. My grandparents have been great rocks which I have clung to. They have always shown courage, and they have always made an effort to make sure I made the right decisions in life. Of my four grandparents, two have passed on, and I miss them very iv much, but I always keep their memories in my heart, and I know they watch over me. My advisors throughout my collegiate career have had drastic influences in my life. Dr. Cremaschi gave me my first work opportunity during my master’s degree, and I learned a great deal from him and my experiences. Dr. Spitler has always been a tough professor, but his lessons are all worth learning, and I have gained invaluable tools throughout my work with him. My graduate advisor, Dr. Fisher, was also my undergraduate advisor. He was the person who opened my eyes to the thermal sciences and created an excitement within me to pursue my master’s degree in this field. His enthusiasm for this field spreads everywhere he goes, and I thank him for motivating me at many times to keep going forward. My master’s degree was not all about hard work, I enjoyed making some great friends throughout this time as well. Sankar, Xiaowei, Steve, Jason, Sam, Michael, Chris and Yang have been there to either play frisbee, soccer, printer tossing, or just to enjoy some dinner with friends. I have always looked forward to their camaraderie. The most influential person in my life is also the love of my life. My wife has been with me for over four years at the time of this writing, and motivated me to continue working no matter how tough it seemed. I have needed her kind words and unwavering love throughout these past years. I have experienced a deep recommitment to Christ since meeting her, and I couldn’t be more proud to call her my wife. v TABLE OF CONTENTS Chapter Page 1 Introduction 1 2 Literature Review 4 2.1 SemiAnalytic Solutions . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.2 Numerical Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.3 Response FactorModels . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.4 ExperimentalWork . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.5 Earth Air Heat Exchanger Analysis . . . . . . . . . . . . . . . . . . . 14 2.6 District Heating . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.7 Soil Property Prediction . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.8 Transport Delay . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.9 Literature Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3 Fluent Reference Model 23 3.1 FluentModel Description . . . . . . . . . . . . . . . . . . . . . . . . 23 3.2 FluentModel Grid Refinement . . . . . . . . . . . . . . . . . . . . . . 24 4 EarthTubeModel 28 4.1 Earth TubeModel Description . . . . . . . . . . . . . . . . . . . . . . 28 4.2 Model Comparisons: VBA Earth Tube vs. Fluent . . . . . . . . . . . 33 4.3 Model Comparison Results: VBA Earth Tube and New Fluent Earth Tube . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 4.4 Earth Tube Conversion to PipeModel . . . . . . . . . . . . . . . . . 40 vi 4.5 Earth TubeModel Guidance . . . . . . . . . . . . . . . . . . . . . . . 45 4.6 Earth TubeModel Conclusions . . . . . . . . . . . . . . . . . . . . . 48 5 MeiModel 49 5.1 Model Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 5.2 Modifications toMeiModel . . . . . . . . . . . . . . . . . . . . . . . 62 5.3 Model Testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 5.3.1 Fully ConstantModel Testing . . . . . . . . . . . . . . . . . . 63 5.3.2 ConstantWall Temperature Testing . . . . . . . . . . . . . . . 63 5.3.3 Initial Condition Independence Testing . . . . . . . . . . . . . 66 5.3.4 High Flow Rate Testing (2DModel Simulation) . . . . . . . . 68 5.4 MeiModel Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . 70 6 PiechowskiModel 73 6.1 Model Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 6.2 Simplified Piechowski Model . . . . . . . . . . . . . . . . . . . . . . . 83 6.3 Analysis of Piechowski Model Grid Dependence Parameters . . . . . . 85 6.4 InitialModel Testing . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 6.4.1 Code Integrity Test . . . . . . . . . . . . . . . . . . . . . . . . 90 6.4.2 Steady State Test . . . . . . . . . . . . . . . . . . . . . . . . . 90 6.4.3 Initial Condition Independence Testing . . . . . . . . . . . . . 94 6.4.4 High Flow Rate Test . . . . . . . . . . . . . . . . . . . . . . . 94 6.5 Comparison to VerificationModel . . . . . . . . . . . . . . . . . . . . 95 6.6 Piechowski Model Conclusions . . . . . . . . . . . . . . . . . . . . . . 102 7 Model Selection 103 7.1 IndividualModel Results Summary . . . . . . . . . . . . . . . . . . . 103 7.1.1 Transfer FunctionMethods . . . . . . . . . . . . . . . . . . . . 103 7.1.2 Earth Tube PipeModel . . . . . . . . . . . . . . . . . . . . . 104 vii 7.1.3 MeiModel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 7.1.4 Piechowski Model . . . . . . . . . . . . . . . . . . . . . . . . . 105 7.2 Selection and Summary of FinalModel . . . . . . . . . . . . . . . . . 106 8 Model Implementation in EnergyPlus 107 8.1 Buried PipeModule Description . . . . . . . . . . . . . . . . . . . . . 107 8.2 EnergyPlus Code Changes . . . . . . . . . . . . . . . . . . . . . . . . 109 8.2.1 WeatherManager . . . . . . . . . . . . . . . . . . . . . . . . . 109 8.2.2 PlantManagers . . . . . . . . . . . . . . . . . . . . . . . . . . 115 8.2.3 Plant Loop Equipment . . . . . . . . . . . . . . . . . . . . . . 115 8.3 Error Handling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 8.4 Preliminary Evaluation of theModel . . . . . . . . . . . . . . . . . . 116 8.4.1 Input File Description . . . . . . . . . . . . . . . . . . . . . . 116 8.4.2 Discussion of Results . . . . . . . . . . . . . . . . . . . . . . . 117 8.4.3 Direct Implementation Comparison: EnergyPlus vs. VBA . . 121 8.5 EnergyPlus Documentation . . . . . . . . . . . . . . . . . . . . . . . 123 8.6 Computation Time . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 9 Conclusions 125 9.1 Model Development . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 9.2 FinalModel Implementation . . . . . . . . . . . . . . . . . . . . . . . 126 9.3 FutureWork . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 A Earth Tube Model Code 138 B Gambit Journal File 143 C Entering Water Temperature Summary 145 D Modified Piechowski Model Code 150 viii E Modified Mei Model Code 168 F Analytic Constant Wall Temperature Code 185 G EnergyPlus IDD Object 187 H EnergyPlus IDF Object 189 I EnergyPlus Input Output Reference Contribution 191 J EnergyPlus Engineering Reference Contribution 195 ix LIST OF TABLES Table Page 2.1 Variable Definitions for equation 2.1 . . . . . . . . . . . . . . . . . . . 6 3.1 Results of Grid Independence Study . . . . . . . . . . . . . . . . . . . 26 4.1 Earth Tube Model Inputs . . . . . . . . . . . . . . . . . . . . . . . . 30 4.2 Variable Descriptions for Eq. 4.1 . . . . . . . . . . . . . . . . . . . . 31 4.3 Variable Descriptions for Eq. 4.7 . . . . . . . . . . . . . . . . . . . . 32 4.4 Calculation of Earth Tube Outlet Temperature . . . . . . . . . . . . 33 4.5 Common Input Data For Earth Tube Models . . . . . . . . . . . . . . 36 4.6 Convection Correlation Information for Circular Pipe Flow . . . . . . 42 4.7 Earth TubeModel Sensitivity Study Results . . . . . . . . . . . . . . 47 4.8 Earth Tube Model Influence Coefficients and Percent Error . . . . . . 48 5.1 MeiModel Update Equation Organization . . . . . . . . . . . . . . . 54 5.2 MeiModel Nomenclature . . . . . . . . . . . . . . . . . . . . . . . . . 55 6.1 Symbols and Variables used in Piechowski Model . . . . . . . . . . . 80 6.2 Piechowski Model Symbol Definitions . . . . . . . . . . . . . . . . . . 83 6.3 Piechowski Input Variable Tabulated Results . . . . . . . . . . . . . . 87 6.4 Simulation Parameters for Piechowski Verification . . . . . . . . . . . 98 6.5 Piechowski Verification: Exiting Temperature RMS Error . . . . . . . 99 6.6 Piechowski Verification: Heat Transfer Rate RMS Error . . . . . . . . 101 8.1 Computation Time Results . . . . . . . . . . . . . . . . . . . . . . . . 124 x LIST OF FIGURES Figure Page 1.1 Three Possible Horizontal Pipe Configurations . . . . . . . . . . . . . 2 2.1 Kusuda and Achenbach Ground Temperature Prediction . . . . . . . 20 3.1 Generic Fluent Buried PipeModel . . . . . . . . . . . . . . . . . . . 25 3.2 FluentModelMesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 4.1 Earth TubeModels: Exit Air Temperatures . . . . . . . . . . . . . . 38 4.2 Earth TubeModels: Exit Air Temperature Deviations . . . . . . . . . 39 4.3 Earth Tube and Fluent Heat Transfer Rates . . . . . . . . . . . . . . 40 4.4 Earth Tube PipeModel Results . . . . . . . . . . . . . . . . . . . . . 43 4.5 Fluent and Earth Tube PipeModel Exit Temperatures . . . . . . . . 44 4.6 Fluent and Earth Tube PipeModel Delta Exit Temperatures . . . . . 44 5.1 MeiModel Domain Description . . . . . . . . . . . . . . . . . . . . . 51 5.2 MeiModel Finite Difference Grid . . . . . . . . . . . . . . . . . . . . 51 5.3 MeiModel Update Equation Locations . . . . . . . . . . . . . . . . . 53 5.4 MeiModel Description: Angular Grid Locations . . . . . . . . . . . . 59 5.5 MeiModel Test Results: ConstantWall Temperature . . . . . . . . . 65 5.6 MeiModel Test Results: Initial Condition Independence . . . . . . . 67 5.7 Mei and Fluent Model Boundary Differences . . . . . . . . . . . . . . 69 5.8 High Flow Test: Mei and Fluent Cross Section Temperature Profiles . 70 6.1 Large Scale View of Piechowski Domain Cutaway . . . . . . . . . . . 74 6.2 Small Scale View of Piechowski Domain Cutaway . . . . . . . . . . . 74 xi 6.3 Case Layout and Node Descriptions for Piechowski Model . . . . . . 77 6.4 Piechowski Radial Grid Distribution . . . . . . . . . . . . . . . . . . 81 6.5 Piechowski Input Variable Study Results . . . . . . . . . . . . . . . . 88 6.6 ConstantWall TemperatureModel Results . . . . . . . . . . . . . . . 93 6.7 Piechowski Initial Condition Independence Results . . . . . . . . . . . 94 6.8 High Flow Rate Test Results: Piechowski Model Cross Section Temperature Distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 6.9 High Flow Rate Test Results: FluentModel Cross Section Temperature Distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 6.10 Piechowski Verification: ExitingWater Temperature Results . . . . . 100 8.1 Main Pipe SimulationManager Routine Flow . . . . . . . . . . . . . 108 8.2 Pipe Initialization Routine Flow . . . . . . . . . . . . . . . . . . . . . 110 8.3 Main Pipe Calculation Routine Flow . . . . . . . . . . . . . . . . . . 111 8.4 Buried Pipe Calculation Routine Flow . . . . . . . . . . . . . . . . . 112 8.5 Buried Pipe Radial Region Routine Flow . . . . . . . . . . . . . . . . 113 8.6 Buried Pipe Water Boundary Routine Flow . . . . . . . . . . . . . . 114 8.7 Results Section 1: 12Feb to 28Apr . . . . . . . . . . . . . . . . . . . 118 8.8 Results Section 2: 31Mar to 8Apr . . . . . . . . . . . . . . . . . . . 119 8.9 Results Section 3: 31Mar to 28Apr . . . . . . . . . . . . . . . . . . 120 8.10 Results Section 4: 2Jun to 17Jun . . . . . . . . . . . . . . . . . . . 120 8.11 Implementation Comparison, EnergyPlus to VBA: Exiting Temperature121 8.12 Implementation Comparison, EnergyPlus to VBA: Total Heat Transfer Rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 xii CHAPTER 1 Introduction Heat transfer from buried horizontal pipes has been studied for a majority of the twentieth century, encompassing such fields as residential and commercial ground source heat pump system design, industrial pipelines, district heating and cooling, buried power cables, and the use of earth tube heat exchangers. These applications have led to studies of soil properties, effects of moisture on underground heat transfer, effects of multiple buried pipes along with different insulation methods. The results of all these studies have shown that the problem is not trivial. There are many factors which increase the difficulty of finding an accurate solution to the problem of the horizontal pipe. The ground is heterogeneous; soil conductivity varies with temperature, soil type, and moisture content. The ground surface changes constantly based on the annual and daily temperature cycles, ground cover materials (grass, concrete slab, etc.), rain, and snow cover conditions. The fluid in the pipe will not maintain constant behavior along the entire length of the pipe, with properties such as viscosity and density changing with temperature. The geometric complexity is increased when multiple pipes are laid in the same trench. Even within this multiple pipe layout there are different configurations, as shown in figure 1.1. These complexities make finding a quick solution nearly impossible. As with any engineering analysis, the solution begins with a set of valid assumptions. The inclusion of a buried pipe model in an annual simulation program is a necessity for properly modeling horizontal ground loop heat exchangers. It is also important for district heating and cooling systems. With lowenergy houses and buildings using 1 Figure 1.1: Three Possible Horizontal Pipe Configurations ground heat transfer, the model becomes even more important, as every degree of precooling required of a water source and every degree of heat lost from a heating pipe can be costly. The main purpose of this current study is to identify a buried pipe heat transfer model which balances accuracy and computation time, and implement that model in the EnergyPlus (2007a) program. The study accomplishes this by the following procedure: 1. Examine previous work in the field. This includes experimental work, which would help define the important aspects of a buried pipe system, and modeling work, which would help define the assumptions that other studies have made in order to form a tractable solution. 2. Develop a reference model. Without experimental data, a robust model must be utilized to provide a comparison tool. 3. Evaluate currently implemented models in EnergyPlus. This step includes looking into models which may not initially be included as buried pipe models, but may be modified to fit such a problem. 4. Develop models for analysis. Each model that seems like a useful option is developed in some form and tested to determine the quality of the results. 2 5. Select a single buried pipe model to implement. The best choice for the model will balance accuracy and computation time, while providing a general purpose solution to fit all applications of buried pipes. 6. Implement the model in EnergyPlus. The model will act as its own module in the EnergyPlus environment, and it will be tested thoroughly to ensure proper operation. 7. Provide user documentation for the model. This documentation will be included as part of the EnergyPlus documentation state what inputs and outputs are available for the model, as well as a reference showing how the model works. During many of the tests and simulations included in this work, hourly weather and entering water temperature data are used as boundary conditions. The weather data are annual hourly data from Atlanta, GA. The entering water temperature data are also hourly data, but it is from an experimental source, the ground source heat pump plant at Oklahoma State. Even though this data is not applicable as model validation, it does provide realistic inputs to the models. This report includes a total of 9 chapters, including this introduction as chapter one. Chapter two is a literature review of different modeling techniques and other aspects of buried pipe design. Chapter three is an introduction to the Fluent reference model used to evaluate the current modeling work. Chapter four evaluates an existing EnergyPlus model, the earth tube model, for modification and use as a buried pipe model. Chapter five evaluates a fully radial model with variable boundary conditions. Chapter six evaluates a model which mixes both cartesian and radial coordinates within the same domain. Chapter seven summarizes results from each model, and selects a model for EnergyPlus implementation. Chapter eight discusses implementation of this model in the EnergyPlus program. Chapter nine summarizes the conclusions of the study and proposes future work. 3 CHAPTER 2 Literature Review In order to find a proper model to implement, it is necessary to look into literature to see what attempts have previously been made. When looking for a buried pipe model, the search may be broadened slightly as well, to include buried ’earth tubes’, since a detailed analysis of a buried air tube shares many aspects with a buried water pipe. The key when looking at these models is to understand the underlying assumptions and ensure that they remain valid for the system at hand. There are several approaches for solving the buried pipe problem, including semianalytic solutions, finite difference solutions, and finite volume solutions. All of these calculation methods are important, but an understanding of experimental work is required in order to validate the applicability of any model. The current search is broad enough to include the experimental studies and understand the important aspects of a buried pipe system design. Some of the other related material comes from district heating literature, which typically includes substantially long piping systems; soil property prediction, which is an important boundary condition for any simulation dealing with the soil; and transport delay, which can be a significant factor, especially with longer piping systems. 2.1 SemiAnalytic Solutions The semianalytic solutions typically consist of a line source or cylinder source solution. In this solution, the soil is assumed infinite in all directions, and an infinitesimally small heat source is applied which represents the pipe. In some cases, this 4 heat source is the input to the model, resulting in pipe outlet temperatures or ground temperatures. The model can be reformulated to instead determine the heat rate when temperatures are given as inputs. Even with this formulation, there are other limitations of the model. The model does not account for ground surface interaction, does not readily account for initial ground temperature gradients, and the line source is also variable along the length of the pipe due to the fluid temperature changing, so this method cannot suitably predict the behavior along the pipe. In many cases these models also use annual average data which can produce a typical annual result. One of the earliest line source models was introduced in a 1948 paper in ASHVE by Ingersoll & Plass (1948). The key in this paper is describing the temperature distribution around a buried pipe. An equation is offered, (2.1), based on the assumptions mentioned above, that can give the temperature of any radial distance away from the line source, valid from the pipe wall to a large distance away. This model assumes that the pipe is the only cause of temperature change in the soil. T − TO = Q 2πk ∞ x e−β2 β dβ = Q 2πk I(x) (2.1) where the variables are defined in table 2.1. The paper offers values of I(x) for numerous values of x, and an equation approximating I(x), which would be more suitable to a simulation program than a lookup table. The original equation was only valid for a small range of x. For the purpose of this current thesis work, the equations were adjusted, and equation (2.2) results, and is valid for a wider range of x. I (x) = 2.303 · log 1 x + x2 2 − x4 8 + 0.1128 · x3 − 0.2886 (2.2) The rest of the paper deals with different pipe situations, and solving examples. The examples range from simple single buried pipes to double buried pipes, and then to simulations where the heat rate changes from month to month. Although these 5 Table 2.1: Variable Definitions for equation 2.1 Variable Definition x √r 2αt T Temperature in soil at any selected distance from the pipe [deg F] To Initial temperature of soil [deg F] Q’ Heat emission of pipe (negative for absorption) [Btu/linear ft/hr] r Distance from centerline of pipe [ft] k Thermal conductivity of soil [Btu/hr/ft2/F] α Thermal diffusivity of soil = k ρCp ρ Density of soil [lb/ft3] Cp Specific heat of soil [Btu/lb/F] t Time since start of operation [hr] β Variable of integration examples do show that the line source approach is a widely applicable method, it still fails in the limitations discussed previously (uniform boundaries, no ground surface interaction). The method still offers incredible simplicity, and has been used for a significant time frame, but with simulations requiring more accuracy, especially in green design, this method does not seem to be the ideal solution. The simplifications for line source applications can be modified and built upon to account for different pipe configurations. In Persson & Claesson (2005), the multipole method is combined with the line source approach in order to account for multiple pipes using a steady state analysis. In this situation, each pipe is considered as a line source, and the resulting temperature distribution around the pipe can be described by an infinite series solution. When this is performed for each pipe, the temperatures and heat flows can be constrained to ’line up’ in between the pipes, thus creating an infinite series of equations describing the system. This infinite series can be truncated 6 to a desired level of accuracy, resulting in a solution accounting for both pipes. In this paper, multiple pipes are simulated in a vertical configuration (one pipe above the other), as well as a horizontal configuration (each pipe next to the other). In the vertical configuration, the supply pipe is simulated both on top and on the bottom. The pipe system is modeled using district heating data (geometry, temperatures) as inputs. It is found that using these simplifications, a fourth order approximation to the infinite solution can provide results of significant accuracy with a computation time of a few seconds. A different analysis using line sources is found in Saastamoinen (2007). In this analysis, a rigorous development of unsteady state equations is performed. The application is a floor heating system, in which there is a conducting layer within the floor where heating pipes are contained. Above this, a layer of carpet is modeled, and underneath the conducting layer is insulation. The development includes a convective boundary condition on the top, and constant temperature boundary condition on the bottom. The unsteady behavior is modeled as a fourier series of periodic onoff cycles. The model assumes the shapes of the isotherms around the line sources are nearly circular near the pipe in order to simplify the model to a tractable solution. The effects of ranging placement of the coils and surface Biot numbers are studied in steady state mode and discussed. Another approach which uses semianalytic methods is performed by Arimilli & Parang (1983) on dual horizontal buried pipes. The method used is a boundary integral method. In this method, a steady state version of Laplace’s equation, represented in Eq (2.3), is written for a twodimensional domain according to boundary integral methods, with the modified version given in Eq (2.4). ∇2T = 0 (2.3) 7 1 2 T (Y ) + B T ∂g ∂n − g ∂T ∂n ds = 0 (2.4) In Eq (2.4), the B subscript refers to integration over the entire boundary, T(Y) represents the temperature at any point Y on the boundary, and g(X,Y) is a harmonic function over the solution domain. The boundary integral involves the temperature along the boundary as well as the derivative at the surface (heat flux for the thermal analysis case). By using this method, the integral can be discretized over a set of segments in the domain. This is similar to a finite difference method, however the boundary integral approach only involves values at the boundaries, and not intermediate values in the center of the domain, which are not of interest in most studies. The resulting equations can be solved with iterative means until convergence. This paper provides the equations required to solve for the general case and describes their simplicity of use, although no validation or verification is given so the accuracy is uncertain. 2.2 Numerical Analysis In many cases, the geometry may be too complex to model properly with a line source approach. For these cases a numerical analysis lends itself well. The numerical approach could be a simple 1D finite difference model, or a 2D finite difference model with variable boundary conditions, or a full 3D finite volume model with unstructured boundary conditions. The approach taken would depend on the pipe configuration, accuracy required, and computation time allowed. The following works span a broad spectrum of numerical solutions and each can provide their own insight into a satisfactory model for implementation in an annual simulation program such as EnergyPlus. There are many variations in ground heat exchangers including concentric tubes, horizontal trenches, and vertical boreholes. In order to develop a model of a vertical heat exchanger, some of the same assumptions can be made, and some different 8 assumptions also. For instance, in a vertical heat exchanger, there will almost always be two tubes in very close proximity. In this configuration, it would be a very bad approach to neglect interaction between the pipes. In many cases, a single horizontal heat exchanger will be placed in a trench, in which case this assumption would be valid. For a vertical heat exchanger, the effects of ground surface will be minimal, as the majority of the pipe will be buried deep underground, while in a horizontal heat exchanger, the entire length of the pipe is affected by the behavior above the ground surface. Chiasson & Spitler (2000) describes the model development of a heated bridge deck using a ground source heat pump system. Lei (1993) presents a validated, computationally efficient model for vertical ground heat exchangers. The model is simplified by neglecting the interface resistance between the pipes and the surrounding material. The effects of the Utube at the end of the pipe are also neglected. The model’s computational efficiency is improved by using a dual domain configuration. A fine region is used in the high gradient regions near the pipe, along with a fully implicit solution scheme. The coarse region farther away is solved using an explicit scheme. The time step used is efficiently selected based on the geometry. In both regions, adaptive grid spacing is employed to avoid wasting computation power where it is not needed. With these assumptions and improvements in place, the cylindrical transient equations are discretized and applied to each node in the domain, and the simulation is solved. The results are compared to experimental results, with high accuracy while the piping system is active. When the flow is stopped, as in the pumping off cycle, the results deviate somewhat. The exiting water temperature during these time steps is not drastically important in most engineering analysis, as long as the model simulates correctly once the flow comes back on. A set of studies is performed by Mei (1988, 1986a) in which both single and double horizontal pipes are simulated. In Mei (1988), a dual pipe is studied, with particular interest in the effects that the second pipe has in the temperature distribution around 9 the first pipe. A quasi3D finite difference code is developed, based on energy balance equations, which simulates both a single pipe and a dual pipe in a tubeovertube configuration. The model assumes that the temperature from inside the pipe up to the pipe wall is radially symmetric, and a circular far field boundary condition is applied which is centered in the bottom tube, and with the top at the ground surface. In order to simplify the model, axial heat transfer is neglected in the pipe and soil, and therefore axial heat is only allowed to flow through the fluid. The model results are compared to experimental data in terms of daily energy absorption in the ground, with the model producing an average error of less than 12%. Mei (1986a) describes three models, each with their own distinct purpose: soil freezing, seasonal ground temperature variations, and thermal interference. Each model is described in detail including mathematical development, computer implementation, validation, and parametric studies. It is found that for soil freezing around the pipe, the effect is minimal as long as the entering temperature is not less than 4◦C, which is typically used as a design constraint anyway to maintain sufficient heat pump capacities. The seasonal ground temperature variation model provides a very high quality simulation of the temperature distribution around a pipe, and also allows for studies of different backfilling materials. This model seems promising for simulation applications and is evaluated further in chapter 5. The thermal interference model is used for any case with multiple buried pipes in the same trench, and is not investigated in the current work, although analysis of multiple buried pipes is recommended as future work. A study which involves both heat and mass transfer was performed by Piechowski (1996, 1998, 1999). In the thesis, Piechowski (1996) described every aspect of the model development in great detail, while the other two papers, Piechowski (1998, 1999), give insight into the model validation and theoretical development, respectively. The model contains aspects that most buried pipe models do not, such as moisture diffusion through the soil, and the addition of thermal storage into the system. The 10 model assumes that the soil temperature a certain distance away from the pipe is not affected by the pipe presence. This becomes the far field boundary condition, and is based on a diurnal and seasonal temperature variation. In order to model moisture transfer in a tractable manner, the mass transfer is only modeled in the vicinity near the pipe, which is where the highest gradients will exist for both heat and mass transfer. The moisture equations therefore have a boundary condition which is in this vicinity. The heat transfer is simulated over a larger domain, with variable boundary conditions. The system is simulated in a 3D finite difference domain, and as with other models, heat and mass transfer are neglected in the axial direction in the pipe material and the soil. This allows for a set of 2D domains to be solved, connected only by the fluid transferred through the pipe. The model is validated against experimental data, and it is found that for a small temperature range, the inclusion of moisture effects (latent heat and moisture migration) does not significantly increase the accuracy, so for most models is would not be necessary. However, the soil thermal conductivity is a very significant input, and improper assumptions of soil properties will lead to inaccurate results. The Piechowski model is developed and tested in chapter 6. Tarnawski & Leong (1993) presented a model of a ground coupled heat exchanger, but with more emphasis on the overall system simulation investigated. The model includes data for several heat pumps, soil properties, heat exchanger configuration (spiral HX included), and pipe materials. The model is verified against simulated results and previously published experimental work, and agrees with a fair amount of accuracy. The verification process was hindered by a lack of information from the verification source, including ground properties and initial conditions. 2.3 Response Factor Models Another heat transfer model type involves the use of transfer functions or response factors. Transfer functions are typically confined to certain systems while response 11 factors can be developed for any system. These are already found in building simulation programs as a way to calculate transient heat transfer through walls using pseudo steadystate means. For each wall construction, transfer functions or response factors are developed which represent how the wall responds if there is a one degree change in temperature on one surface. Then superposition can be applied to determine the effects of general temperature changes. There are also response factors used in ground heat transfer, but there are limitations to this method. A short time step response factor model for vertical borehole heat exchangers was developed by Yavuzturk & Spitler (1999) and validated by Yavuzturk & Spitler (2001). In this model, a series of response factors are developed using a transient two dimensional finite difference system. With the response factors calculated, the model is able to produce the temperature response of the borehole field to a step heat input from the ground heat exchanger. This model type is appealing because once the response factors have been generated, the responses are computationally fast. The response factors developed in this work are also valid for short time steps, which is suitable for an energy simulation program such as EnergyPlus. This model proves to work quite well for vertical heat exchangers. Horizontal heat exchangers have a disadvantage for such an approach in that they are greatly affected by the ground surface. For a vertical heat exchanger, only the small section at the top will be affected, so the response factors are allowed to ignore this behavior. For a shallow horizontal heat exchanger, the response factors would need to include the response in the ground to not only a step heat input from the heat exchanger, but also to a step change in temperature or heat from the ground surface. The response factors must also be able to handle modeling a certain depth of soil. Currently, there is a horizontal surface heat exchanger model available in the EnergyPlus program that is based on a modified model from Strand (1995). This model uses transfer functions which are modified to include a heat source within the construction. This is also found when modeling radiant floor systems, in which 12 case a warm pipe is embedded within the construction. This method is quite useful and again, it is very fast compared to iterative finite difference techniques (the transfer functions in place are also iterative due to required surface heat balances). The major problem with using this model as a general purpose buried pipe model is the depth. Existing transfer functions do not work with very heavy constructions, and with approximately 1m of soil, the EnergyPlus transfer functions diverge. Response factors can be developed for any system such as a horizontal heat exchanger, however the work is beyond the scope of this thesis. The development of response factors is recommended as future work. 2.4 Experimental Work To ensure that any model is relevant and useful, there must be some experimental work in the field. This experimental work can provide a means of validating a model, provide ranges of useful input data, characterize important properties of the system that must be investigated and accounted for in the model, and also provide a proof of concept to ensure that model development is worth the time and effort. Inalli & Esen (2004) provided details on a horizontal ground source heat exchanger system in heating season in Turkey. The study provides key data which could be used as input to a new model development. The study also provides experimental evidence that the soil properties play a significant role in the design of the heat exchanger, as it can easily lead to oversizing of the components and therefore the whole system. In this study, the average coefficients of performance (COPs) of the system (2.66 2.81) were lower than in other experiments due to oversizing of the equipment. The data provided is presented more in terms of seasonal and daily averages. The data provided is not suitable as a source for model validation. An article that discusses not only experimental apparatus, but also useful design techniques, was presented by Klimkowski et al. (1985). The paper studies the long 13 term effects of average Uvalue for a given horizontal heat exchanger. Data is measured, and then resolved back into average conductance. After a period of continuous operation of the heat pump, the Uvalue approaches a steady constant, one value in winter, and another value in summer. It is recommended to use this value as a design criteria. The paper then offers a design process, which simply utilizes average earth temperature, fluid temperature, and a given Uvalue, along with a required capacity, to calculate a required HX length in both summer and winter. The longer of these two results is selected as the required pipe length. This type of analysis is a very simplified model which may provide a decent prediction of total annual energy usage, but gives no insight into the daily or hourly operation. These smaller time steps, along with accounting for energy storage in the ground, is where savings may be found using a more robust model, especially when looking at low energy applications. 2.5 Earth Air Heat Exchanger Analysis When developing a model for a buried pipe, it is desirable to look to different aspects of the heat transfer field in order to possibly find a suitable model that can be modified to fit the needs of a buried pipe. One example is an earth tube. In an earth tube, air is passed through a buried duct in order to precool the air in summer, and preheat the air in winter. In this way, the demands of the conditioning equipment are reduced, while the ventilation requirements are easily met with fresh air. The earth tube models differ in the assumptions, as air is passing through the pipe instead of water or glycol. However, if a model shows good performance under earth tube conditions, it may be useful to modify it and analyze it further as a buried pipe. AlAjmi et al. (2006) presented an earth tube model which employs a variable soil temperature variation, and applies steady state behavior and an effectiveness correlation in order to calculate exiting temperature conditions and overall heat transfer. The model is verified against previously published data and proves a high level of 14 accuracy. The model is then implemented into the TRNSYS environment and used in a full building simulation to show the benefits of using an earth tube, as the cooling load is decreased substantially. Baxter (1992, 1994) provided a full description of experimental results of an earth tube system in both heating and cooling modes. The results show that one of the main factors that affects the earth tube performance is the temperature gradient that forms in the axial direction near the pipe. This indicates that assuming constant soil temperature near the pipe should be used cautiously, as it could be erroneous. Bojic et al. (1999) simulated a pvc and a steel pipe together, and showed results of the model. The model was based on a 3D Finite Element Method, although no details were given about the meshing specifics. While the results are given, no validation or verification is offered. One conclusion is that the earth tube operates better in the summer months than the winter months. This may be due to the simulation climate and other temperatures used in the simulation. Another approach to modeling an earth tube is given by De Paepe & Janssens (2003). In this paper, the authors developed a thermo hydraulic design of earthair heat exchanger. In this model, both the heat transfer effects and the fluid mechanics effects are modeled, resulting in a very detailed study. The purpose of this model is to evaluate the balance between pressure drop and heat transfer when designing a heat exchanger. The paper develops an NTU correlation for the earth tube heat exchanger. A ’specific pressure drop’ is developed to allow a designer to then calculate required lengths of parallel tubes in order to achieve the required capacities. A very simplified model of an earth tube, which is implemented in the EnergyPlus simulation program, is developed by Lee & Strand (2006). This model is developed in a very simplified manner in order to achieve two main goals. The first goal is implementation in an annual simulation program. A full finite volume study of an earth tube would be too computationally intensive to recalculate every time step 15 of an annual simulation. The second goal is available input data. A simulation engine such as EnergyPlus requires that the models have input to which designers have access. The model takes in simple average annual data, and at each time step, calculates a soil temperature in the near pipe region, and calculates soil, pipe, and convective resistances. A steady state approach then uses an energy balance between the soil and air to create an equation to calculate exiting air temperature as a function of pipe geometry, fluid characteristics, and soil data. The model can also account for fan power if mechanical ventilation is selected. This model is an almost instantaneous calculation model, and therefore lends itself well to being converted to a buried pipe model and reanalyzed. This model is evaluated in detail in chapter 4. Mihalakakou et al. (1994) provides a model which discusses both heat and mass transfer around the pipe, and accurately predicts the exiting air temperature. The model uses a 3D finite element approach which uses concentric rings as the elements. The model offers the conservation equations which are applied to each node and include both heat transfer and moisture transfer. With boundary conditions specified, the model was used to simulate an earth tube over a period of 15 days. Over this time, experimental data was measured, and the results of the model were in very satisfactory agreement. This model was implemented within the TRNSYS environment, and provides a means of properly predicting earth tube exiting temperature as well as soil temperature variations. In a study by Tzaferis et al. (1992), two sets of experimental data were measured, and eight different earth air heat exchanger models were simulated under these experimental conditions, and their results were compared, along with the model sensitivity. The models were split in two categories: those which analyzed ground heat transfer, and those which only analyzed heat transfer through the pipe. In the latter of these, temperatures at or near the pipe wall were required as inputs. The models range from assuming constant pipe wall temperature (analytic integral) to using finite difference 16 schemes. The root mean square (RMS) error for nearly all the models were below 3.5%, showing that even simplified models can provide decent results for an earth tube. 2.6 District Heating One of the larger scale applications of a buried pipe model is in district heating systems. In these systems, a long pipe loop connects a large number of buildings, providing hot water, which can then be tapped and used as a primary heat source or as a source for heat pumps. In these systems, a good model for the pipe heat transfer effects should be used since even a small percentage of energy savings can be quite significant. A simplified model by Bohm (2000) takes an approach quite different from other models. The model assumes that at each time step, there is a soil temperature that can be used in steady state equations to mimic transient behavior. This temperature would account for the storage effects in the ground by using a temperature from a previous time step. In this manner, the soil temperature near the pipe would be ’lagged’ by a certain amount. The position of this undisturbed temperature in the ground is found by numerical simulations and experimentation. A function is developed for this temperature, although several restrictions were placed on the analysis. For the verification, a constant pipe temperature was used, and it is recommended to only maintain constant model inputs for very brief periods. In this study, as with most district heating systems, the pipe configuration is a dual pipe, with the pipes laid side by side. The results do compare well with experimental data taken, although, some of the model inputs come from predicted, or experimentally determined data. A model of an indirect district heating system is offered by Chuanshan (1997). This model looks deeply at very detailed effects in a district heating system, and not so deeply at predicting the heat transfer of the buried piping. The study does offer 17 details on model parameters which are very useful if a buried pipe model is developed and tested in a district heating simulation. Eriksson & Sunden (1998) detailed some new aspects of modeling, while still accounting for both heat and mass transfer in the soil. The model appears to be a 1D approach, with the solution using a tridiagonal matrix algorithm (TDMA) approach. The key to this particular model is the investigation of pipe material breaking down through time. The thickness and conductivity of the pipe changing over time is studied over a longterm simulation, and the results show that the thermal conductivity of the pipe increases over time, but does not make a significant impact in the overall heat exchanger performance. 2.7 Soil Property Prediction In order to properly simulate a ground coupled heat exchanger, the soil effects must be modeled accurately as well. The effects of surface temperature variation is even more important when dealing with a horizontal heat exchanger than with a vertical heat exchanger because the pipe or duct is near the surface for the full length, not just a small initial section. This section discusses models and information found for predicting not only ground temperatures, but also ground properties. A study by AbuHamdeh & Reeder (2000) shows that the thermal conductivity is highly dependent on other soil properties such as density, but the application of the study is not relevant for a buried pipe calculation. Fischer (1983) gives a brief summary of fourteen soil heat and moisture transfer models, describing each model’s capabilities and limitations, but does not include any model development itself. A model in Kusuda & Achenbach (1965) is one of the more popular models for predicting ground temperatures. The model is based on an assumption that the earth is a semiinfinite solid, and that the surface boundary condition can be described by a sinusoidal temperature variation. In this manner, a simple correlation is given, based 18 on some simple input data, but the ground surface temperature and average annual amplitude of temperature variation must be known. The main equation (2.5) uses an exponential term to account for depth into the ground (delay), and a cosine term to account for annual variation. T = A − BO · e −x· √ π DT cos 2πθ T − x π DT − PO (2.5) In equation (2.5), T is the ground temperature at a given time θ and depth x. D is the thermal diffusivity of the soil, T is the time period of the annual variation, and should be the same units as θ. PO would be a phase angle to account for the date of minimum surface temperature. A is the annual average surface temperature, and BO is the annual average variation in surface temperature. As an example, to show the effects of ground heat diffusion, figure 2.1 shows the annual ground temperature for a series of depths. In this example, A was set to 18.333◦C, BO was set to 11.667◦C, D was set to 1.765E3 m2/hr, PO was set to zero, and since the run was an annual simulation and diffusivity was in hours, the constant T was set to 8760 hours. Kumar & Kaushik (1995) discusses a soil temperature model which is based on Kusuda & Achenbach (1965), but includes soil properties (conductivity, specific heat, density, and diffusivity) for ten different soil consistencies (dry light, damp heavy, average rock, etc.). 2.8 Transport Delay One aspect that is typically not dealt with in detail is the transport delay problem. Most simulation programs have pipe models, but do not account for the time required for the fluid to pass through the pipe. In the past, this may have been neglected due to the additional computational strain that this would place on a simulation program. With computers becoming faster, this aspect should not be ignored when it is applicable to include it. 19 Figure 2.1: Kusuda and Achenbach Ground Temperature Prediction 20 Hanby et al. (2002) develops a model that is based on second order equations used at several nodes along the length of a pipe. The model accounts for the thermal capacitance of the fluid and the conduit wall. Each node is a wellmixed zone, and the node separation distance is based on capacity values and average fluid residence time. This residence time is a factor which can be assumed based on flow behavior or precalculated using CFD studies. The model shows very close agreement with experimental data. Saman & Mahdi (1996) details some previous models that attempted to predict the transport delay problem, along with where these models fell short. The model in the paper is based on analytic energy balance equations discretized along the length of the pipe, and solved using finite difference methods. The results do follow the expected behavior, but no validation is offered in the paper, so the accuracy of the model cannot be determined. 2.9 Literature Summary This chapter has been a review of literature involving buried pipe heat transfer models, soil property prediction, and transport delay models. The buried pipe models included semianalytic methods, numerical simulation methods, and response factor methods. The response factor methods reviewed were generated for vertical heat exchangers, which allowed the ground surface to be ignored. The generation of response factors for a shallow horizontal heat exchanger is not in the scope of this thesis work, but is highly recommended as future work. The semianalytic methods required a large number of assumptions and typically required input such as the line source strength; doing so requires a guess and therefore introduces high uncertainty into the results. The numerical methods varied in computational requirements and accuracy. The models selected for development in this work are the Earth Tube model, the V. C. Mei model, and the Piechowski model. These three 21 give varying degrees of assumptions and accuracy, and should give a nice comparison of what work has been performed, and which model is worthy of development within the EnergyPlus system. The soil property prediction by Kusuda has been used for years covering different fields of study. For the models currently developed in this work, the Kusuda and Achenbach method will be used to predict boundary condition temperatures and surface temperatures when appropriate. The transport delay models will currently not be developed, as the main goal of this thesis work is to get a buried pipe heat transfer model into EnergyPlus. The introduction of a transport delay model could be obtained by linking several smaller pipe models together once one is developed. This is left as another good recommendation for future work. However, with a transient buried pipe model, the effects of a transport delay are somewhat captured. The basic pipe model allows the inlet temperature to pass directly through to the pipe outlet. With the addition of any transient effects, the immediate transport will encounter a delay. Further investigation will reveal if the models need to be developed together to create one dualpurpose model. In conclusion, three numerical heat transfer models will be developed which use the Kusuda and Achenbach soil temperature prediction method, and do not include an actual transport delay aspect. These three models will be compared against a verification model discussed in chapter 3, and a decision will be made to implement one of the models weighing accuracy vs. computation time. 22 CHAPTER 3 Fluent Reference Model 3.1 Fluent Model Description When creating or testing any model, there must be a method of ensuring the model is operating properly. In some cases, analytic solutions may be available which can be used as a validation source. Another way of performing validation is through experimental results. If one can create an experimental apparatus or facility to which the model can be applied, the results can be compared to prove the model’s validity. In the current case, experimental data is not available, and so direct validation is not possible. In some limiting cases, analytic means are used in this work to prove model operation. These cases require substantial assumptions to allow an analytic solution to be developed. For the majority of this current work, verification will be performed. In verification, the model is compared against a more robust, or previously validated model. Fluent CFD software is employed for the current situation. In most cases, Fluent is used to simulate fluid flow in diverse situations. However, the finite volume engine within Fluent is powerful enough to simply switch off the momentum equations, and only solve the conservation of energy equation in solid regions. This allows conduction to be modeled easily within Fluent, along with convection flow in a tube, all in the same model. Another benefit of using Fluent is the flexibility of boundary conditions. There are features within Fluent that allow users to input their own functions that can be spatially or temporally dependent, and applied to any boundary of the model. The functions are written in the C programming language, can be easily hooked to a 23 model, and can access many internal variables to Fluent’s solver. This is very useful when simulating, for instance, the annual variation of ground surface temperature. The user defined function can reach within Fluent, retrieve the current simulation time step, and return a value of surface temperature for that particular hour of the year. For some situations where the boundary condition is spatially dependent, such as the Kusuda and Achenbach boundary (a function of time and depth), the function can reach into Fluent and get the depth coordinate in order to return the appropriate earth temperature for that time and depth. The other boundary conditions can all be applied as needed, some symmetric, some known boundary temperature. The exact configuration will depend on the model being verified. Figure 3.1 shows a typical Fluent model for this work. The buried pipe channel along with the relevant boundary conditions, are shown. The boundary conditions will vary depending on the current model being tested. Throughout this work, a number of simplified models will be compared to Fluent models that have been carefully generated to closely match the simplified model. The details of these comparisons are discussed in each individual model chapter. 3.2 Fluent Model Grid Refinement In order to develop a valid Fluent model, an independent grid must be demonstrated. By using the gambit journal file, it is easy to reproduce models quickly. An example of the journal files used in this study can be found in appendix B. In order to ensure that the grid is independent, the grid was reproduced several times, each with increasing refinement of the grid. The grid was refined in both the pipe and soil regions, each with its own given node spacing value. For the grid independence check, a steady state situation was utilized, with simple boundary conditions. The ground surface was held at 310K, while the other spatial boundaries were adiabatic. The pipe inlet was air at a uniform 1 m/s and 260K. The entire domain was initialized to the conditions 24 Figure 3.1: Generic Fluent Buried Pipe Model 25 Table 3.1: Results of Grid Independence Study at the pipe inlet, and the simulation ran to convergence. Convergence was kept at the default residual values (1E03 for most variables, 1E06 for energy equation). After convergence was attained and the results were checked for satisfactory convergence (no oscillation, etc.), the total heat transfer rate on the pipe wall was calculated by Fluent. This was used as the grid independence monitor. The grid was simulated for four different mesh conditions, and table 3.1 shows the results of the runs. Table 3.1 shows that the results are in fact converging on a solution as the mesh becomes more and more refined. The table also shows that the time required for the simulation to run is becoming much larger. In order to provide quality results in a reasonable time frame, a compromise must be made between accuracy and computation time. For further computations, the mesh for case 3 (10 minute run time) is selected. This simulation gives quite good results with only about 5% improvement over the previous run. The mesh is shown in figure 3.2. The grid independence is one of the most important checks to make when performing finite difference studies. Other studies to verify the Fluent model under different circumstances will be performed on a modelbymodel basis throughout this work. 26 Figure 3.2: Fluent Model Mesh 27 CHAPTER 4 Earth Tube Model 4.1 Earth Tube Model Description One of the models which is used to predict ground heat transfer from a buried pipe was developed for an Earth Tube system (Lee & Strand 2006). In this earth tube system, ventilation air is brought into an underground duct from a certain distance away from the zone or conditioning system, and the air is tempered, either preheated or precooled depending on the season, and then brought into the zone or through the conditioning equipment. This pretreatment of the air helps to reduce the overall heating and cooling equipment loads while satisfying fresh air requirements. The model makes many assumptions, but requires minimal input data and has a very short computation time. These features lend themselves well to first stage design work in an annual simulation program. Because of the many assumptions, this model seems to be one of the simplest models to apply to a buried pipe. As with many models, one must begin by understanding the underlying assumptions, in order to ensure that the model is valid for a particular system, and not oversimplified. • The first assumption is that fluid flow within the pipe is thermally and hydrodynamically developed. This assumption can be easily made when the length of the pipe is much greater than the diameter. Developing flow will only occur within a short initial section of the pipe, and will not greatly influence the thermal behavior when the pipe is long. This assumption allows a single correlation 28 to be used for any given fluid and flow regime. • The next assumption is that soil temperature in the pipe vicinity can be calculated using the model developed by Kusuda & Achenbach (1965). This earth tube model makes use of the soil temperature model to calculate a soil temperature one annulus from the center of the pipe. • The temperature in the soil near the pipe is not affected by the pipe itself. This assumption allows the soil temperature to be uniform along the entire axial length of the pipe. This allows the solution of the system to be greatly simplified, but may impact the accuracy. • The soil is homogeneous and maintains constant thermal properties (thermal conductivity) throughout any time step, although at each time step, the properties can be adjusted to a new temperature dependent value and held constant yet again. • The fluid is axisymmetric at any cross section, and the pipe maintains a constant cross sectional area in the axial direction. • At any given time step, the heat transfer behavior can be modeled as steady state. In order to be useful, the model inputs must be available. Typically, this requires simplifying the model to accommodate data that most people can access. This particular model requires input of properties such as system geometry, material properties, fluid flow properties (including a convection correlation), and annual surface temperature data. The full set of input data is described in table 4.1. The system geometry should be available to a designer, and the fluid flow properties can be approximated based on design conditions. Soil properties and annual surface temperature data must be estimated. The soil properties will be variable based on moisture content, 29 Table 4.1: Earth Tube Model Inputs Pipe Data Geometry: Buried depth, Axial length, Radius, Wall thickness Properties: Pipe wall thermal conductivity Soil Data Properties: Thermal diffusivity, Thermal conductivity Other: Annual average surface temperature Annual surface temperature amplitude Fluid Data Flow: Flow rate, Convection correlation Properties: Density, Specific heat, Prandtl number Viscosity, Thermal conductivity and the surface temperature data will be location specific. Soil properties may be determined from an onsite thermal conductivity test, and surface temperature data can be approximated from annual weather data. Once the model has the required input data, it uses hourly outdoor air temperature and a series of equations to calculate the air temperature exiting the earth tube. These equations are described here. The first equation calculates the ground temperature in the vicinity of the buried pipe. This equation is developed from the work by Kusuda & Achenbach (1965), but arranged to fit the attributes of this particular model: Tgnd,z = Ts − As ∗ e −z∗ √ π 365∗α ∗ cos 2 ∗ π 365 ∗ t − tshift − z 2 ∗ 365 πα (4.1) The variables used in equation (4.1) are defined in table 4.2. 30 Table 4.2: Variable Descriptions for Eq. 4.1 Symbol Description Tgnd,z Ground temperature at pipe depth [C] Ts Annual average surface temperature [C] As Annual surface temperature amplitude [C] z Pipe burial depth [m] α Ground thermal diffusivity [m2/day] t Simulation Run Time [day] tshift Day of minimum surface temperature (i.e. 0 = 1Jan) [day] With the ground boundary temperature found at the pipe depth, focus can be switched to the fluid. To define the fluid to pipe interface, a convection coefficient is required. The correlation (not listed here) is a function of the working fluid (air, water) and the flow regime (laminar, turbulent). For typical fluids such as air and water, there are many correlations available in common heat transfer and fluid flow textbooks. Whichever correlation is used, the outcome is the convection coefficient, h, which can then be used in further calculations. The soil data, pipe data, and fluid data have now all been initialized and steps can be taken to compute the exiting fluid temperature. First, the thermal resistances of each material and interface are calculated, then the total resistance, and therefore the total conductance, can be calculated with the following equations. Convection resistance : Rc = 1 2πr1h (4.2) Pipe wall resistance : Rp = 1 2πkp ∗ ln r1 + r2 r1 (4.3) Soil resistance : Rs = 1 2πks ∗ ln r1 + r2 + r3 r1 + r2 (4.4) Total resistance : Rt = Rc + Rp + Rs (4.5) Total conductance : Ut = 1 Rt (4.6) 31 The paper then applies a heat balance to a differential length of pipe, where the heat transfer to the air is equated to the heat storage within the air. Because of the simplifications that previously resulted in a uniform temperature along the length of the pipe, the outlet air temperature is a simple expression which is a function of conductance, fluid properties, entering fluid temperature, ground temperature near the pipe, and pipe length. For simplification, the paper defines an intermediate term that can be calculated, which simplifies the final exiting temperature expression, and is given as equation (4.7). A = m˙fCf ∗ ln Tf,i − Tgnd,z − Ut ∗ L m˙fCf (4.7) The variables used in equation (4.7) are described in table 4.3. Table 4.3: Variable Descriptions for Eq. 4.7 Symbol Description m˙f Mass flow rate of fluid Cf Specific heat capacity of fluid Tf,i Entering fluid temperature Tgnd,z Ground temperature near pipe L Pipe length Using this intermediate term from equation (4.7), the outlet temperature can be calculated based on different scenarios regarding entering fluid temperature, Tf,i, and ground temperature, Tgnd,z, as described in table 4.4. As this description shows, this model is a very simple model to implement, and lends itself well to annual simulation programs, where it may be called thousands 32 Table 4.4: Calculation of Earth Tube Outlet Temperature Situation Equation Tf,i > Tz,t Tf,ex (L) = Tz,t + eA Tf,i < Tz,t Tf,ex (L) = Tz,t − eA Tf,i = Tz,t Tf,ex (L) = Tz,t of times during one simulation. When implementing this model in a spreadsheet environment, the annual simulation at an hourly time step took less than one second to complete. The simplicity does come at a price, however, and that price is accuracy. The assumption that the ground temperature near the pipe is not affected by the presence of the pipe goes against what most ground source heat pump developers count on, and that is heat storage in the ground. This model will be tested against more rigorous models to determine how much accuracy is lost in the assumptions and simplifications made. 4.2 Model Comparisons: VBA Earth Tube vs. Fluent The earth tube model, as implemented in VBA, requires the inputs found in table 4.1. To reproduce the model in Fluent, modifications to the Fluent model discussed in chapter 3 are required. Geometrical constraints such as pipe diameter, depth, length, and thickness are all very similar. In the VBA implementation, there is no finite difference solution, so the domain is simply the pipe along with a representative soil temperature with which the pipe interacts. For the Fluent model, however, the geometric constraints must include distances to apply boundary conditions in the 3D domain. These distances represent the width and height of the soil that is modeled. The soil is modeled underneath the pipe, and around the side of the pipe. The VBA earth tube model does not take into account any heat transfer from these regions. In an effort to produce the most direct comparison of the two models, the most 33 appropriate way to implement these boundary conditions in the Fluent model is to use adiabatic boundary conditions on these faces. The soil and fluid properties may be directly placed within both models, which helps ensure agreement between the two. The fluid to pipe convection coefficient will cause some difference between the models. For the earth tube model, a heat transfer correlation is used which is based on Reynolds number, and therefore based on bulk flow rate through the pipe. In Fluent, the local velocity and boundary layer will be simulated, and so a uniform convection coefficient will not be modeled. The last item that changes is the soil temperature prediction. The earth tube model uses annual average surface temperature data to produce a single correlation as a function of depth and time. This then produces a temperature in the soil near the pipe at each time step of the simulation. In Fluent, this would over specify the problem, since the soil is being modeled as a solid volume, not just the near pipe vicinity. The inputs to the earth tube model are the annual average surface temperature and the annual surface temperature amplitude. These two inputs can also be used in Fluent to develop a userdefined transient boundary condition. Some details on userdefined functions are presented in section 3.1. The surface temperature can be approximated as a sine wave throughout the year with the same given average and amplitude. In summary, the models can be made very similar, but with three significant differences: 1. The soil is modeled from the surface to a distance below the pipe and in the region to the side of the pipe in Fluent. The ’deep ground’ and ’farfield’ boundaries will initially be modeled as adiabatic to account for the fact that the VBA earth tube model ignores heat transfer interaction with these regions. 2. The air flow through the tube is modeled as a finite element mesh in order to predict heat transfer behavior from the air in Fluent. In the earth tube model, a correlation is used based on bulk flow through the tube. 34 3. The soil surface boundary condition will be a userdefined sine wave approximation of a yearly cycle with a given average and amplitude in Fluent. In the earth tube model, the same average and amplitude are used in a correlation to predict the temperature of the soil in the nearpipe vicinity, which then becomes the boundary condition for the model. In other words, the earth tube calculations only use the surface temperature to precalculate the soil temperature at the pipe depth for any given day of the year, and the surface is not directly modeled. Note that the VBA implementation was not directly verified against the EnergyPlus implementation. A paper by Lee & Strand (2006) covers the development of the model and implementation of the model in EnergyPlus. The work in this paper was used as the guideline for the current VBA model implementation. One major difference that may occur between the VBA and EnergyPlus implementation is the time step management; the VBA implementation uses a constant one hour time step, while EnergyPlus varies the time step as needed throughout the simulation. The model is updated with the current ground temperature and the earth tube is simulated. Based on a verification of the governing model equations, the model presented in the paper and the VBA model should produce the same results. 4.3 Model Comparison Results: VBA Earth Tube and New Fluent Earth Tube With the proper mesh developed for the Fluent model, the earth tube model as implemented in VBA and the Fluent model can be directly compared. The detailed differences between the models can be found in Section 4.2. The inputs used which are common to both models are given in Table 4.5. The inputs which are not shared between the models are summarized below. For each input, the method used by each model is given. 35 Table 4.5: Common Input Data For Earth Tube Models Input Description Value Units Pipe Data Burial Depth 2.5 [m] Length 25.0 [m] Radius 0.25 [m] Thickness 0.01 [m] Thermal Conductivity 200 W m−K Soil Data Thermal Conductivity 1.5 W m−K Fluid Data Working Fluid Air [−] Flow rate 0.785 m3 s Simulation Data Surf Temp Avg 15 [C] Surf Temp Ampl 5.6 [C] • Ground Surface Temperature – VBA Earth Tube: Use average surface temperature and amplitude data along with equation (4.1) to generate soil temperature near pipe. – Fluent Earth Tube: Use average surface temperature and amplitude data as a transient sine wave: Tsurf (t) = Tsurf,avg + Tsurf,amplsin (A ∗ t). The A coefficient is modified to fit the time unit conventions in Fluent. This equation is similar to the full Kusuda & Achenbach (1965) correlation, but does not require the exponential because the depth is zero at the ground surface. 36 • Deep Ground Heat Transfer – VBA Earth Tube: This model does not take into account any heat transfer interaction to deep ground or heat transfer to the side of the pipe. – Fluent Earth Tube: Soil surrounding pipe is modeled, but all surfaces are set as adiabatic to mimic current earth tube behavior. • Fluid Properties – VBA Earth Tube: This model uses constant values of air properties such as density and specific heat. – Fluent Earth Tube: Fluent contains libraries of material properties, some of which are constant, and some of which are calculated by equations of state (such as ideal gas). – Note: For this model comparison, the air will not see drastic temperature or property changes, so this change is not a significant source of error. The main output of the VBA earth tube model is the outlet temperature from the tube. The earth tube model inherently assumes steady state behavior. During an annual simulation, the ground temperature and inlet temperature are updated, and a new steady state solution is found. For both models, a time step of one hour is used. The simulation is then operated for 8760 time steps to include a full annual simulation. In the earth tube model, hourly weather data is read directly. For the Fluent model, the entering air temperature is a boundary file, which consists of ’ifthen’ structures to select the outdoor dry bulb temperature for any given hour of the year. The weather data used is from Atlanta, GA. Although Fluent can be used to output any data desired, the variable of most interest is the exiting fluid temperature. This is what will be compared to the earth tube model. The goal of this comparison is to determine how well the earth tube 37 Figure 4.1: Earth Tube Models: Exit Air Temperatures can predict ground heat transfer, compared to the more robust and thorough Fluent model. Other comparisons are performed as needed. The exiting air temperatures are plotted for the entire simulation year in Figure 4.1. The numerous simplifications inherent in the earth tube model result in some disagreement with the Fluent model. The deviations between the models are plotted in figure 4.2. The earth tube model is a steady state model at each time step, while Fluent solves the transient equations of momentum and energy for each time step. At the beginning of the simulation, Fluent initializes the entire domain to a steady condition based on the values at the initial time step. For the entire year, the trends of both models agree, but the actual values do not. For most of the year, the models deviate by about 5◦C. In the summer, the Fluent model shows better precooling capabilities, while in the winter, it shows better preheating capabilities. The earth tube model appears to be in error, but in a conservative manner. This seems to be due to the better ground heat storage modeling in Fluent. Another way to characterize the error between the two models is to show the different heat transfer rates predicted. This is similar to presenting the pipe delta 38 Figure 4.2: Earth Tube Models: Exit Air Temperature Deviations temperature (inlet to output), but gives a more applicable aspect to the results. Most designers will care more about the heat transfer rate than predicted temperatures. Figure 4.3 shows the heat transfer rates predicted by the two models. The heat transfer rates between the two models show the same basic trends, but differ drastically in actual values. Through the simulation, the heat transfer rate predicted by Fluent is up to four times as high as the earth tube model prediction. The model predicts higher heat transfer rates in both heating and cooling conditions. This shows that the earth tube model is, in fact, a conservative model. When a designer is using the earth tube model in EnergyPlus within a simulation, it should be noted that the earth tube may actually perform better than the model predicts. Although the earth tube model does deviate from the Fluent model, the extreme increase in simplicity for model implementation and computation speed still make it a possible candidate to become a buried water pipe model. A Fluent model will also be developed to verify the model accuracy when the earth tube is converted to a buried pipe model. This model development is described in section 4.4. 39 Figure 4.3: Earth Tube and Fluent Heat Transfer Rates 4.4 Earth Tube Conversion to Pipe Model This section describes the changes made to the earth tube model in order to simulate a buried pipe with water flow. For simplicity in this study, many simulation parameters will remain constant. The surface temperature will still be modeled by a sinusoidal wave in Fluent, using the same average and amplitude as the earth tube model. The pipe depth will remain constant, although the pipe diameter will change. The new pipe will be a 0.1 [m] nominal diameter pipe. The other boundary conditions will remain constant except for the inlet conditions. For this, a new set of boundary values based on ground source heat pump water temperatures will be used. This data was recorded at the Oklahoma State University Heat Pump Laboratory. The data is the entering water temperature to a set of vertical boreholes, so the data cannot be used as validation for a horizontal heat exchanger loop. However, the data does provide a semirealistic input to the model. There is one key issue with simulating a short section of pipe and not an entire loop. The error in exiting temperature will have a tendency to propogate through the entire loop, and the entering water tem 40 perature will be a function of the exiting water temperature. By using this entering temperature as the boundary condition, the small amount of error in the loop is corrected each new hour. If the full loop were run along with a heat pump model, this error may accumulate. This simulation is out of the context of the Fluent reference model, and requires long term validation data to account for these long term effects. However, as the current model is only for a buried pipe, using this entering water temperature data will ensure that the simulation is performed in the most realistic manner possible without drastic changes to the model itself. The values are summarized in appendix C, and the same values will be used in the Fluent model and the earth tube model. The water flow rate will also be different from the air flow rate. The water velocity will be decreased from 1.0 m/s to 0.15 m/s. For a pipe of 0.1 m diameter, this is on the order of 0.0012 m3/s, or 19 GPM. For typical standard conditions, this results in a Reynolds number of approximately 12,000. For the water flow, Fluent will automatically solve the momentum and boundary layer equations using water data instead of air data once this is selected as the working fluid. For the earth tube model, a new convection correlation will be required. With this Reynolds number known, the convection correlation can be specified in equation (4.8), with the the correlation applicable over the range specified in Table 4.6, as found in Incropera & DeWitt (2002). NuD = 0.027Re4/5 D Pr1/3 μ μs 0.14 (4.8) For the case of this simulation, the minimal property variations will not have a drastic effect on the results because the temperatures throughout the entire simulation vary only within a small range. This assumption should allow the last term in equation (4.8), μ μs 0.14 , to be neglected, thus simplifying the simulation. With this convection correlation in place, the simulation can be performed. For an annual simulation with the given data, the resulting temperatures are shown in figure 41 Table 4.6: Convection Correlation Information for Circular Pipe Flow Conditions: Turbulent Flow Fully Developed 0.7 ≤Pr≤ 16, 700 ReD ≥ 10, 000 L/D≥10 Variable Description ReD Reynolds number based on diameter Pr Fluid Prandtl number μ Viscosity at fluid bulk temperature μs Viscosity at film temperature 4.4. There are three curves shown: Pipe Inlet Temp, Pipe Outlet Temp, and Temp Deviation. The inlet and outlet are basic flow temperatures, although on this plot, the differentials are actually hard to see. With this simulation, the temperature difference from inlet to outlet is small, on the order of 1◦C. This is why the differentials are also plotted. The temperature difference is plotted on a different scale to make it easier to visualize. As can be seen on the plot, the temperatures fluctuate widely. This is simulating the heat pump demand during the week, and then off on the weekends. The model predicts about 1◦C of precooling in the summer, and 0.5◦C of preheating in the winter. The Fluent model was also modified, as discussed previously. The main change in the mesh involves a smaller pipe diameter. This smaller pipe requires a tighter mesh as well, but this is a simple task to perform once a baseline journal file has been created for Gambit. The new pipe mesh was placed in Fluent, and the same simulation that was performed for the earth tube was repeated, with two changes: Pipe inlet temperature and working fluid. The simulation was initialized to data 42 Figure 4.4: Earth Tube Pipe Model Results at the first time step until steady state, and then the transient annual analysis was performed. The exit temperatures for both the Fluent pipe model and the earth tube pipe model are shown in figure 4.5. Figure 4.5 shows that the Fluent model predicts a much higher precooling rate during the summer, and a significant but lower preheating rate during the winter. To better visualize the differentials, the absolute value of the deviations is plotted in figure 4.6. This difference means that, as with the earth tube model comparison, the earth tube is a conservative simulation. If used as a design tool, it appears that the actual performance will be higher. This is most likely due again to the Fluent model capability to better model heat storage in the ground. For the earth tube, the air has an effect on ground temperatures near the pipe, but this effect is small relative to the effect from a water pipe. For this water pipe, the heat storage in the ground is more prominent. Improper modeling of the ground capacity along with previously discussed model assumptions led to poor results from the modified earth tube model. 43 Figure 4.5: Fluent and Earth Tube Pipe Model Exit Temperatures Figure 4.6: Fluent and Earth Tube Pipe Model Delta Exit Temperatures 44 4.5 Earth Tube Model Guidance The results of the pipe study show that the earth tube model does not give good results when using water as the working fluid. One reason for this is the added thermal storage in the system. In a typical earth tube, the air flows through quickly, and the effects of the air are not as prominent in the ground as a flow of water. The model creates a steady state simulation at each time step, assuming a constant temperature in the ground near the pipe. With water, the ground absorbs and rejects significant heat which causes the variation in temperature along the length of pipe to be more intensified. An improvement to this model could be found in the calculation of the ground temperature in the near pipe region. Currently, surface temperature data is used in a correlation and this results in a temperature at the pipe depth, which is used as a boundary condition to the model. This temperature, however, is too close to be considered a farfield boundary. A serial model is not robust enough, and there should be some feedback from the pipe when generating a ground temperature. This is a difficult improvement when the model remains steady state. The earth tube model is built upon deep underlying assumptions. Although the model has limitations, some general guidance is provided for situations where the earth tube model should provide higher accuracy. Due to the constant temperature assumption along the length of the pipe, the earth tube model accuracy should be improved in designs where the earth tube will not encounter major temperature variations through the ground along the length of the tube. This is similar to the statement that the earth tube model will be more accurate for cases where ground temperatures are dominated by groundsurface interaction rather than groundpipe interaction. • Minimal earth tube pipe length: The shorter the tube, the less chance there is for a temperature variation within the tube, and therefore within the ground. 45 When a shorter tube is in place, however, the entrance effects are much more definitive. The heat transfer coefficient will be higher in this entrance region due to the developing boundary layer, but overall, the heat transfer will be reduced. • Increased flow rate: When the flow rate is reduced, the temperature difference in the fluid will be greater. This axial variation poses an opportunity for inaccuracy in the model. Therefore, it may be inferred that the model will increase accuracy with higher flow rate. • Soil properties: With the earth tube model assuming steady state behavior at every time step, higher accuracy should be found in simulations that better approximate steady behavior. Higher values for soil conductivity and diffusivity lead to faster heat diffusion through the soil. This will in turn reduce the possibility for temperature variation in the soil. This will become a case that performs similar to a steady situation. With regard to soil properties, the higher diffusivity used should result in the higher accuracy model. In addition to the model input guidance given above, a parametric study was performed with the earth tube model within the EnergyPlus simulation program to determine which of the model inputs have a significant impact to the simulation output. For the inputs with higher sensitivity, more time should be spent increasing the quality. If an input shows negligible sensitivity, time will be wasted trying to improve its quality. A simple earth tube model was set up for the annual simulation. The simulation contains a full building simulation, and the earth tube is applied to a single zone. The representative output for the model is selected as the furnace heat transfer rate. One at a time, the variables shown in table 4.7 were adjusted, and the change in coil energy requirements were recorded. The sensitivity was nondimensionalized according to equation (4.9) with inputs X and outputs Y. The subscripts refer to 46 either 0 for the base case simulation or i for the simulation of case i. Sensitivity = Yi−Y0 Y0 Xi−X0 X0 (4.9) With equation 4.9, the study revealed the data in table 4.7. Table 4.7: Earth Tube Model Sensitivity Study Results Case Variable Baseline New Inp Input Chg Output Chg Sensitivity 0 Baseline      1 Flow (m3/s) 1 1.25 0.25 0.0007 0.0028 2 Radius (m) 0.3048 0.4048 0.3281 0.0003 0.0008 3 Length (m) 50 75 0.5 0.0063 0.0127 4 Depth (m) 3 4.5 0.5 0.0002 0.0003 In addition to testing the sensitivity of the model, influence coefficients can provide another level of meaning to how the model responds to different inputs. Application of influence coefficients to building simulation are discussed in Spitler et al. (1989). These coefficients are similar to the sensitivity discussed previously. However, the estimated percent error is developed which includes the estimated error (uncertainty) in the input value. Equation (4.10) is used to define the type 2 dimensional influence coefficients, and equation (4.11) is used to define the percent error. Influence Coefficient = Y0−Y1 Y0 X0 − X1 (4.10) Percent Error = Influence Coefficient · Estimated Error in Input · 100 (4.11) The influence coefficients defined here are dimensional as they contain the inverse units of the input parameter. When the percent error is calculated, the estimated 47 error in input carries the input units also, which results in a dimensionless percent error, as expected. The four case results are presented in Table 4.8 in terms of influence coefficients and percent error. The uncertainty of the inputs was set to 5%. Table 4.8: Earth Tube Model Influence Coefficients and Percent Error Case Variable Influence Coefficient Percent Error 1 Flow (m3/s) 0.00276 1.378% 2 Radius (m) 0.00251 1.253% 3 Length (m) 0.00025 0.127% 4 Depth (m) 0.00012 0.058% This influence coefficient study shows that the percent error in results is highest in the flow rate and pipe radius. The sensitivity of the model was shown above to be higher in the flow rate and pipe length input parameters. With this information, it is assumed that the highest possibility for model improvement is in the flow rate estimation. The simulation indicates that the earth tube model is sensitive to the flow rate, and a small uncertainty in the input parameter produces a percent error of about 1.4%. 4.6 Earth Tube Model Conclusions The results of this chapter indicate that the earth tube model, as it operates currently, is not a suitable candidate for the buried pipe model. Other models will be compared throughout this study, most of which include a more substantial set of input data and require significantly more computation time. If none of these seem suitable, the earth tube model may become the basis of a new model that better handles heat storage in the soil and pipe. 48 CHAPTER 5 Mei Model 5.1 Model Description Mei (1986a) developed a model to predict heat transfer around a buried horizontal pipe. This model uses an explicit finite difference scheme, and develops update equations, which are used in a pseudo3D environment along the pipe. As one would expect, this model uses symmetry to split the soil region in half, centered along the pipe. This model is considered simple due to a single cylindrical coordinate system. This gives a less complex algorithm within the code compared to other models which may mix coordinate systems. The original code, however, seems more complex, as the model was programmed using older FORTRAN conventions. This model was translated to VBA and tested in the Excel environment. There are actually numerous advantages when using this model. One, already mentioned, is the single coordinate system. This simplifies the heat transfer equations in the domain. Another advantage this model has is that it allows for a backfill material to be used around the buried pipe. This ability is already in the code, so the only items a user must input are the material properties of the backfill material, and how thick the material is, and the program will automatically simulate the region accordingly. In this same manner, if no backfill material is used, this region can simply be set to the outer soil properties, and the program will run as normal, with soil filling the backfill region. There are also some disadvantages of this model. The single radial coordinate system is centered on the centerline of the pipe, and extends the radius to the ground 49 surface (although the domain can be fully underground). The model layout can be seen in figure 5.1. This layout indicates that there will be a single node at the ground surface at each cross section. This does not seem to be as effective as spreading out nodes along the ground surface. However, testing may indicate that a single node is suitable for representation. Another disadvantage may be the distribution of grid points. This model does allow the grid points to be nonuniform in the radial direction in the soil, with fixed spacing in each of the backfill and pipe wall regions. However the variation in grid spacing can lead to stability problems. The radial coordinate system is applied to nodes in different ways throughout the domain. The water at any given pipe cross section is represented as a single node. The pipe wall is represented by radial nodes; variations in the angular (θ) direction are ignored. This correlates to the assumption that the heat transfer up to the outer pipe wall is axisymmetric. Therefore, at a given cross section and radius, the pipe temperature is assumed constant for all angles. Moving outward radially, the pipe outer wall interfaces with the backfill region. The code has an interesting effect through this backfill region. There are nodes set up as if a full 2D (radial and angular) grid was set up, but the heat transfer only passes in the radial direction. This is analogous to breaking the domain up into a 3D grid, but ignoring the heat transfer along the axial direction. This allows a firststage analysis of heat transfer varying with angle and therefore depth. Moving outward the next section is the main soil section. In this region, the domain is fully 2D, both angular and radial heat transfer effects are modeled. To better illustrate this grid, figure 5.1 describes the various domains and figure 5.2 shows the layout of the finite difference grid. As mentioned previously, there are a number of assumptions that are used by the model, in order to develop a tractable solution. Some of these may turn out to be invalid in certain situations. The main assumptions used, as presented by Mei, are summarized as follows: 50 Figure 5.1: Mei Model Domain Description Figure 5.2: Mei Model Finite Difference Grid 51 • Homogeneous soil, which allows constant properties within the soil region. However, with the addition of the backfill region, a secondary homogeneous material (either a second soil type or pipe insulation) can be supplied for the analysis. • Fluid temperature and velocity are constant at any cross section. This allows a single node to represent the fluid temperature and velocity within the pipe at any axial point. • Pipe depth is great enough that the domain boundary can be considered a farfield boundary. This assumption allows the model to ignore surface heat transfer effects directly. • There is only one pipe in the ground in the given domain. • Heat transfer within the pipe up to the pipe wall is axisymmetric. This allows the heat transfer within this region to be 1D radial. The model operates on a dual time step. The pipe wall and water nodes run at a shorter time step in an attempt to ensure convergence, while the backfill and soil nodes run at a much longer time step. This allows the accuracy and stability of a short time step only in the regions where it is required. In the backfill and soil regions, where the temperature gradients are smaller, a longer time step is more suitable. As one would expect, since the model is broken down into four different regions (water, pipe wall, backfill, and soil), there are four different temperature arrays which hold the data. The water temperature array is defined in the program as U (Length, Time Step). The pipe wall temperature array is defined as T (Radius, Length, Time). In the backfill region, the heat transfer is modeled as pseudo 2D; the temperature can vary angularly, but the heat transfer along the angular direction is ignored. The array will therefore have four variables, and is defined as F (Radius, Length, Angle, Time). The final temperature array is the soil temperature 52 Figure 5.3: Mei Model Update Equation Locations array. It is defined as S (Radius, Length, Angle, Time). These arrays are updated sequentially, according to specific update equations for the simulation domain and model conditions. The flow of the model in its entirety is similar to other models. The model reads parameters supplied by the user, initializes the domain, starts the transient iterations, updates boundary temperatures, updates the water and pipe wall temperatures on the smaller time step and then updates the soil and backfill regions on the larger time step, writes output at each larger time step, and repeats until the simulation is complete. More detail is given later in this section as required. The update equations differ based on the region and conditions. For instance, the pipe wall to water interface equations change depending on whether the flow is on or off. The different locations to which an update equation might be applied are shown in figure 5.3. Table 5.1 references the locations shown in figure 5.3, describes the conditions, and references the equation number of the update equation discussed in the following paragraphs. In the equations, the independent variables have been simplified. The time step is not listed as a subscript as are the other variables. Instead, the convention has been 53 Table 5.1: Mei Model Update Equation Organization Location Number Description (Figure 5.3) (in text) (Conditions) A 5.1 Water Node, Flow ON, J=1 A 5.4 Water Node, Flow OFF, All Other J B 5.3 Inner Pipe Wall, Flow ON B 5.5 Inner Pipe Wall, Flow OFF C 5.6 Internal Pipe Wall Node, All Conditions D 5.7 Outer Pipe Wall, All Conditions D 5.8 Inner Backfill Node, All Conditions E 5.9 Internal Backfill, All Conditions F 5.10 Outer Backfill, θ = 0 F 5.11 Outer Backfill, 0 < θ < π F 5.12 Outer Backfill, θ = π G 5.13 Soil, θ = 0 G 5.14 Soil, 0 < θ < π G 5.15 Soil, θ = π H 4.1 Farfield, All Conditions 54 Table 5.2: Mei Model Nomenclature Parameters and Subscripts Region Extents Pipe Length Inlet to Outlet J 1 to NZ Angular Location Top to Bottom of Domain M 1 to NRAD Pipe Wall Radial Location Pipe Centerline to Outer Pipe Wall I 1 to NPIPE Backfill Radial Location Outer Pipe Wall to Outer Backfill I 1 to NFZ Soil Radial Location Outer Backfill to Outer Domain Boundary I 1 to NSOIL Example Variable Description UJ Water Temp at Length J TI,J Pipe Wall Temp at Radius I, Length J FI,J,M Backfill Temp at Rad I, Length J, Angle M SI,J,M Soil Temp at Radius I, Length J, Angle M applied which uses the prime ’ notation to denote the updated time step. Aside from the time step, the independent variables are given in table 5.2. As mentioned, the variables in this section represent both the updated temperatures and the temperature from the previous time step. The updated temperatures are denoted with a prime, such as F I,J,M. Note that with this dual time step program, the main outer time step, which runs the soil and backfill regions, is denoted Δt1. The inner time step, which runs the coil wall and water nodes, is denoted Δt2. These variables will appear throughout these equations. In this domain, there are a number of radii which are very specific. The general 55 radius term is denoted rI , and can be used within any of the regions. There are also four specific radii that represent specific interfaces in the domain. These, along with the last variables of interest in these equations, are denoted and defined in the following: • r1: Inner Pipe Wall Radius • r2: Outer Pipe Wall Radius / Inner Backfill Radius • r3: Outer Backfill Radius / Inner Soil Radius • r4: Outer Soil Radius / Farfield Boundary Radius • h: Water to Pipe Convection Coefficient • V : Flow Velocity • Δr: Radial spacing: This variable has a constant value in the pipe wall and backfill regions, but is nonuniform in the soil region. Therefore, the radial spacing is dependent on the current radius. The first set of equations presented here are applied to the water and inner pipe wall nodes when the flow is on. U 1 = 1 − V dT dz − 2hΔt2 ρfCpf r1 U1 + V dT dz E.W.T. + 2hΔt2 ρfCpf r1 T1,1 (5.1) U J = 1 − VdT dz − 2hΔt2 ρfCpf r1 UJ + V dT dz UJ1 + 2hΔt2 ρfCpf r1 T1,J (5.2) T 1,J = 1 − αpipeΔt2 Δrpipe 2 Δrpipe + 1 r1 − 2αpipeΔt2h Δrpipekpipe T1,J + αpipeΔt2 Δrpipe 2 Δrpipe + 1 r1 T2,J + 2αpipeΔt2h Δrpipekpipe UJ (5.3) 56 Equations (5.4) and (5.5) are used when the flow in the pipe is off. U J = UJ + 2kpipe (T2,J − T1,J)Δt2 ρfCpf r2 1log 1 + Δrpipe r1 (5.4) T 1,J = U J (5.5) Once the simulation has left the inner pipe wall surface, the equations are no longer dependent on the flow being on or off. Therefore, all of the remaining equations in this discussion are applied for both flow on and flow off conditions. Equations (5.6) and (5.7) are for the rest of the pipe wall nodes, both interior nodes and outer pipe wall surface node, respectively. In equation (5.7), a term TAVG has been introduced. Within the pipe wall, there is a single set of pure radial nodes. Outside the pipe wall, there are multiple nodes in the angular direction, even though angular interaction is neglected. Therefore at the outer pipe wall interface, the grid turns from a single 1D radial grid to a number of radial grid points. The single interior node must interact with the multiple grid points outside the pipe wall. Therefore, it is assumed that the one outer pipe wall node interacts with the average of all innermost backfill nodes. This average of the backfill nodes which are directly adjacent to the pipe wall is called TAVG for brevity. T I,J = 1 − αpipeΔt2 Δrpipe 1 Δrpipe + 1 rI − αpipeΔt2 Δr2 pipe TI,J + αpipeΔt2 Δrpipe 1 Δrpipe + 1 rI TI+1,J + αpipeΔt2 Δr2 pipe TI−1,J (5.6) T Npipe,J = 1 − 2 kfill kpipe αpipeΔt2 ΔrpipeΔrfill − αpipeΔt2 Δrpipe 2 Δrpipe − 1 r2 × TNpipe,J + 2 kfill kpipe αpipeΔt2 ΔrpipeΔrfill TAVG + αpipeΔt2 Δrpipe 2 Δrpipe − 1 r2 TNpipe1,J (5.7) 57 When the simulation has reached this point, both the water and the pipe wall nodes have been simulated. This occurs at a time step of Δt2. The simulation then moves on to the backfill region and simulates at a time step of Δt1. There are three subregions in the backfill region, the inner node, interior nodes, and the outer node. Each of these are defined individually in equations (5.8) through (5.12). F 1,J,M = T J,NPipe (5.8) F I,J,M = 1 − αfillΔt1 Δrfill 1 Δrfill + 1 rI − αfillΔt1 Δrfill 2 FI,J,M + αfillΔt1 Δrfill 1 Δrfill + 1 rI FI+1,J,M + αfillΔt1 Δrfill 2 FI−1,J,M (5.9) Once the simulation has reached the outermost backfill node, the grid has three possible equations based on the node location. Even though the update equations in the interior of the backfill region are only derived in the radial direction, the soil region is derived in both the radial and the angular directions. The grid points along the interface of the backfill and soil regions are derived to match the soil region. Because of the angular dependence, the grid points at the angular domain extents will encounter the adiabatic (symmetric) boundary. This leads to three different update equations, one for M=1 (θ = 0 radians at the top of the domain), one for M=NRAD (θ = π radians at the bottom of the domain), and one for the interior points. These locations are described in figure 5.4. The equations follow. 58 Figure 5.4: Mei Model Description: Angular Grid Locations F NFZ,J,1 = 2ksoilαfillΔt1 kfillΔrsoil,1Δrfill S2,J,1 + 2αfillΔt1 Δrfill 1 Δrfill − 1 2r3 FNFZ−1,J,1 + αfillΔt1 (r3Δθ)2 2FNFZ,J,M+1 + FNFZ,J,1 × 1 − 2ksoilαfillΔt1 kfillΔrsoil,1Δrfill − 2αfillΔt1 Δrfill 1 Δrfill − 1 2r3 − 2 αfillΔt1 (r3Δθ)2 (5.10) 59 F NFZ,J,M = 2ksoilαfillΔt1 kfillΔrsoil,1Δrfill S2,J,M + 2αfillΔt1 Δrfill 1 Δrfill − 1 2r3 FNFZ−1,J,M + αfillΔt1 (r3Δθ)2 (FNFZ,J,M+1 + FNFZ,J,M−1) + FNFZ,J,M × 1 − 2ksoilαfillΔt1 kfillΔrsoil,1Δrfill − 2αfillΔt1 Δrfill 1 Δrfill − 1 2r3 − 2 αfillΔt1 (r3Δθ)2 (5.11) F NFZ,J,NRAD = 2ksoilαfillΔt1 kfillΔrsoil,1Δrfill S2,J,M + 2αfillΔt1 Δrfill 1 Δrfill − 1 2r3 FNFZ−1,J,M + αfillΔt1 (r3Δθ)2 2FNFZ,J,M−1 + FNFZ,J,NRAD × 1 − 2ksoilαfillΔt1 kfillΔrsoil,1Δrfill − 2αfillΔt1 Δrfill 1 Δrfill − 1 2r3 − 2 αfillΔt1 (r3Δθ)2 (5.12) With these equations applied at M=1 through M=NRAD, the backfill region outer interface is complete, and the simulation can move to the soil region. The soil domain is meshed into a full 2D polar grid. This mesh is similar to the interface of the backfill and the soil in that there will be three different update equations based on the angular location. These three equations account for the adiabatic surfaces at the top and bottom of the left side of the domain, as well as the internal points. This mesh structure can be visualized as before in figure 5.4. Equations (5.13), (5.14), and (5.15) are for M=1 points, Internal points, and M=NRAD points, respectively. 60 S I,J,1 = 1 − 2αsoilΔt1 Δrsoil,IΔrsoil,I1 − αsoilΔt1 Δrsoil,I1 2 Δrsoil,I1 − 1 rI − αsoilΔt1 (rIΔθ)2 × SI,J,1 + 2αsoilΔt1 Δrsoil,IΔrsoil,I1 SI+1,J,1 + αsoilΔt1 Δrsoil,I1 2 Δrsoil,I1 − 1 rI SI−1,J,1 + αsoilΔt1 (rIΔθ)2 2SI,J,2 (5.13) S I,J,M = 1 − 2αsoilΔt1 Δrsoil,IΔrsoil,I1 − αsoilΔt1 Δrsoil,I1 2 Δrsoil,I1 − 1 rI − αsoilΔt1 (rIΔθ)2 × SI,J,M + 2αsoilΔt1 Δrsoil,IΔrsoil,I1 SI+1,J,M + αsoilΔt1 Δrsoil,I1 2 Δrsoil,I1 − 1 rI SI−1,J,M + αsoilΔt1 (rIΔθ)2 (SI,J,M+1 + SI,J,M−1) (5.14) S I,J,NRAD = 1 − 2αsoilΔt1 Δrsoil,IΔrsoil,I1 − αsoilΔt1 Δrsoil,I1 2 Δrsoil,I1 − 1 rI − αsoilΔt1 (rIΔθ)2 × SI,J,NRAD + 2αsoilΔt1 Δrsoil,IΔrsoil,I1 SI+1,J,NRAD + αsoilΔt1 Δrsoil,I1 2 Δrsoil,I1 − 1 rI SI−1,J,NRAD + αsoilΔt1 (rIΔθ)2 2SI,J,NRAD−1 (5.15) Thus all interior points of the calculation domain have been acknowledged. The only points which do not have an update equation given here are along the outer farfield boundary. The depth of each boundary point in this radial domain is calculated 61 based on input geometry. With the current node depth and simulation time known, the Kusuda and Achenbach correlation which has been discussed previously in this work is used to update the node at each outer time step. The Kusuda correlation used here is very similar to equation (4.1). More detail is provided in section 5.3.4. 5.2 Modifications to Mei Model After the initial implementation of the Mei model, some modifications were required to fit better with the current implementation language and some features were removed or modified. The structure of the original program was very difficult to understand because the programming did not include some common structures such as the ’Else If’ command. Without a structure such as this, the code must wind in and out with hundreds of ’GOTO’ statements, rendering it quite difficult to follow. Thus the program was modified and updated with newer command structures whenever possible. As for actual model features, there were two items which were modified: specification of the backfill region, and temperature arrays to store previous values of backfill and soil nodal temperatures. In the original program, the use of a backfill region was optional. This allowed a user to place the pipe in direct contact with the soil. However, this option introduces many additional equations. In order to simplify the program, the option was removed, and the backfill region is now in place all the time. This does not actually limit the program in any way. To simulate the pipe without a backfill region, one can simply specify the backfill properties to be the same as the soil properties. The second change involved the storage of values from previous time steps. Although it may have been a simple omission, the code in Mei’s report never included storing the backfill and soil temperatures at the end of each time step. This was revealed during testing as discussed in section 5.3.3. The addition of storage resulted in a model that performed as expected. 62 5.3 Model Testing After initially debugging the model, it appeared to be operating correctly. However, tests were performed to help reveal any underlying model problems. The first test was a simple test where the entering water temperature and all boundary conditions were set to the same temperature. If the model were operating correctly, there would be essentially no response from the model, just constant temperatures everywhere. For the second test, the inputs to the Mei model were adjusted to mimic a constant pipe wall temperature condition. This was set by carefully adjusting the conductivities and other parameters. In the third test, the model was run for different initial conditions to determine convergence characteristics and to ensure that the model was properly updating each temperature node. In the fourth test, the flow rate was increased to eliminate axial heat transfer effects and to assess the 2D effects calculated by the models. 5.3.1 Fully Constant Model Testing This test was designed to ensure that the model algorithm and update equations are defined properly. If the update equations are not implemented properly, heat transfer may be simulated through nodes when it should not, simply because of an error in a coefficient. During this test, the boundary temperatures and inlet water temperature were set to a constant value of 50◦F. The model responded properly to these boundary conditions. This test served to establish some confidence in the model implementation. 5.3.2 Constant Wall Temperature Testing This test approximates an analytic solution. The analytic solution is simply a constant pipe wall temperature solution for exiting fluid temperature as a function of entering water temperature. A more detailed discussion of the analytic solution is 63 given in section 6.4.2. In this test, the Mei model was modified to approximate a constant pipe wall. The conductivities in the backfill region, soil region, and pipe wall were increased to the order of 999 Btu/hrftF, and the outer boundary temperature was held constant. With the high conductivities, this boundary temperature should transfer quickly to the inner pipe wall. This test was performed with a boundary temperature of 50◦F and an entering water temperature of 0◦F. The test was repeated for several different pipe lengths to see the model behavior. At each run, the pipe inner surface temperature along the length of the pipe was recorded. With this numerical scheme, it is nearly impossible to get the pipe temperatures to exactly 50◦ while achieving model convergence. Therefore, an average surface temperature for each run was found, and this was used as the surface temperature in the analytic solution. The exiting water temperature for the numerical case was then compared to the analytic solution. The main analytic equation governing the exiting fluid temperature is given here as equation (5.16), but as mentioned, the full development is in section 6.4.2. As is also mentioned in that section, the analytic solution was applied to small segments of the pipe in order to improve the comparison with the Mei model. T2 = k0.023 ρDV μ 0.8 Pr0.4πL Ts − T1 2 + ˙ mCpT1 k 20.023 ρDV μ 0.8 Pr0.4πL+ ˙ mCp (5.16) The results of this study for the Mei model are given in figure 5.5. At the lengths reported, the models agree to within 5%. This test removes the soil modeling from the simulation, and with this high agreement, increases the confidence in the water modeling of the simulation. The dual time step, which allows the soil to run at a long time step while the pipe wall and water node can run at a short time step, is required for convergence. With a thin pipe, the pipe node spacing becomes very small. When these small nodes are 64 Figure 5.5: Mei Model Test Results: Constant Wall Temperature 65 encountered, a substantially short time step must also be used. During this specific test of constant wall temperature, conductivity values were set above realistic values. This was required in order to quickly translate the boundary temperatures to the pipe wall. The conductivity plays a role in the update equation coefficients, and when these values were increased, the model had serious convergence problems. In order to achieve convergence, the model was run at a 30 minute outer time step, and a .001 minute inner time step. For normal conductivities, initial tests revealed that for an outer time step of one hour, an inner time step of 0.1 minute was required for convergence. The minimum system time step in EnergyPlus is one minute. In order for this model to be implemented in EnergyPlus, an outer structure would be required to simulate this model a number of times during a single EnergyPlus time step. This stability limit does not completely eliminate the possibility of implementing the Mei model in EnergyPlus, but it does introduce another level of complexity to the model and the implementation. 5.3.3 Initial Condition Independence Testing In order to ensure that the model was properly updating the temperatures at each node, a series of initial conditions were placed on the model. The model was then run for a period until the results converged. Figure 5.6 shows the results. In the first test, an initial condition of 0◦F was placed on every node in the model domain. In the second test, the entire model domain was initialized to the Kusuda and Achenbach predicted temperature at the pipe burial depth. In the third test, the model domain was initialized to 100◦F. Figure 5.6 shows that even with the extreme initializations of tests one and three, the model converges, though over two weeks worth of simulation time was required! 66 Figure 5.6: Mei Model Test Results: Initial Condition Independence 67 5.3.4 High Flow Rate Testing (2D Model Simulation) The Mei model does not account for the axial heat transfer effects in the soil or pipe wall. The water in the pipe is the only region which carries heat along the length of the pipe. The Fluent model, however, includes heat transfer in all dimensions. Because of this difference, it is necessary to perform a test which can simulate the ability of the models to simulate in just two dimensions (eliminating the axial dependence). This is performed by substantially increasing the flow rate through the pipe. With this higher flow rate, the fluid will experience almost no temperature change along the entire length of the pipe. If the boundary conditions also do not change along the length of the pipe, two dimensional heat transfer will result. This test marks the first comparison between the Mei model and the Fluent verification model. The Mei model is a cylindrical domain. The Fluent verification model is a rectangular based domain. For many of the boundary conditions, this will not have an effect. In the farfield boundary condition, however, a deviation will occur between the two models. Both models will use the Kusuda and Achenbach correlation to develop the boundary temperatures. In order to create a pseudo steadystate model, the simulation time will be set to zero so that this Kusuda temperature field will not change throughout the iterations. The Fluent model will use the Kusuda and Achenbach directly on the ground surface (no outdoor air interaction), and the farfield (side wall) boundary, and the deep ground boundary. The Mei model will use the Kusuda and Achenbach correlation directly on the cylindrical perimeter of the model. The Mei model does not differentiate between the deep ground and side boundaries, because it is one outer surface boundary in cylindrical coordinates. Figure 5.7 shows the boundaries for both models. While increasing the flow rate, a value will be reached where the axial effects are 68 Figure 5.7: Mei and Fluent Model Boundary Differences negligible. In the Mei model, the fluid enters at 37◦C, with a flow rate of 0.24 m3 s . The flow was increased similarly in the Fluent model. The removal of axial heat transfer from the model should help to reveal two things: whether the models agree when only 2D heat transfer effects are concerned, and if the different boundary type seem to have an effect on overall model agreement. Both of these models are transient in nature, and therefore include the effects of heat storage in the soil and other regions where appropriate. In order to perform this analysis, the models were run until a steady condition was achieved. In the Fluent model, one could turn the model to a steady behavior by switching model parameters, and letting it run until convergence. However, in the interest of comparing these two models in their base configurations, they were both left as transient. To monitor for steady conditions, some variables had to be recorded. For the Mei model, representative temperatures at numerous locations throughout the domain were printed at each time step. In the Fluent model, temperature output is not as convenient, so the overall heat transfer from the pipe wall was recorded, and monitored for steady 69 behavior. Figure 5.8: High Flow Test: Mei and Fluent Cross Section Temperature Profiles Figure 5.8 shows the temperature distributions for both models. As the figure shows, the temperature distributions do not agree well. The temperature distribution along the farfield boundary follows a behavior that would be expected from the Kusuda and Achenbach correlation, but the rest of the domain shows deviations. The Mei model has already shown problems with initial condition independence and stability. The inaccuracy of the 2D calculation is another serious limitation of the Mei model. 5.4 Mei Model Conclusions The Mei model was introduced as a model with intermediate complexity. The single cylindrical coordinate system indicated a model which could produce quality results while balancing the computation time that may be encountered with a more complex model. However, throughout the initial testing, the model did not produce qual 70 ity results and encountered several problems which eliminate it as a candidate for implementation within EnergyPlus: • Stability • Accuracy • Ground Surface Interaction • Initial Condition Independence The first problem, stability, is a result of two things: an explicitly formulated model, and widely varying grid spacing. This varying grid spacing is partially accounted for by the dual time step method used. However, in the initial testing there were several simulations that were unusable because instability caused divergence. This is the main problem of the Mei model. The second problem involves the accuracy of the model. This first set of initial tests, where the domain is idealized in one form or another, show that the Mei model does not agree well with the verification model. The results are expected to worsen as the domain is placed back to its default state where the boundary conditions are both spatially and temporally dependent. The third problem is that the model does not directly interact with the ground surface. The farfield boundary condition of the model is simply a correlation of depth and time. The problem with this is that the model never takes into account the short time scale behavior of the outdoor air. With the other problems encountered during testing, a solution for this problem was not investigated. The fourth problem is that the model showed a very slow independence from initial conditions. The model took almost fourteen days to ’forget’ the initialization of the model. A ’warmup’ period would be required for the model in addition to the base warmup periods that EnergyPlus uses. These problems indicate that theMei model is not currently the ideal model for implementation into EnergyPlus. In chapter 6, the Piechowski model will be evaluated, 71 and depending on these results, chapter 7 will reveal which model is implemented into EnergyPlus. Also, if no model is suitable in its current state, a decision will be made regarding major modifications to one of the models in order to improve the quality, so that it can then be implemented. 72 CHAPTER 6 Piechowski Model 6.1 Model Description Piechowski (1996) developed a quasi3D finite difference buried pipe heat transfer model. Although a set of 3D grid points is generated, the effects of heat transfer from one section of pipe length to the next is only modeled within the fluid. In other words, the soil and pipe material do not model axial heat transfer. Piechowski assumes that the larger heat transfer gradients occur in a radial direction away from the pipe surface, with relatively smaller gradients in the axial direction. Therefore, at each pipe cross section, a 2D grid is developed. The grid is based on a cartesian coordinate system, with boundaries at the ground surface, a deep ground depth, and along the farfield surface on the side. The domain splits the pipe along the middle of the cross section in order to accommodate symmetry. Within this cartesian grid, there is one singular node at each cross section that includes the soil near the pipe, the pipe, and the fluid within the pipe. This layout is shown in figure 6.1. Using a single node alone to represent the pipe would not provide sufficient accuracy, and so a secondary coordinate system was developed, centered upon this representative node. This secondary domain is in radial coordinates. In this radial system, where the grid spacing is more refined than in the outer cartesian system, both heat and mass transfer are calculated, while in the larger grid, only sensible heat transfer is accounted for. The smaller grid layout is shown in figure 6.2. Within the program, there are three distinct temperature arrays. Each of these accounts for a different section of the simulation domain. The first temperature array 73 Figure 6.1: Large Scale View of Piechowski Domain Cutaway Figure 6.2: Small Scale View of Piechowski Domain Cutaway 74 is T (TimeStep,Pipe Length,Soil Depth,Soil Width). This temperature array has a value for the entire Cartesian soil domain. In this region, the grid is a 3D Cartesian coordinate system, and includes a representative value for the rectangular node which contains the pipe. This node will be discussed further in a later section. The second temperature array is held in the variable TR (TimeStep,Pipe Length,Radius). This temperature array represents the values within the radial nearpipe region. In this section, the heat transfer is assumed to be axisymmetric, and the grid is therefore 1D (radial) throughout each cross section. When this region is analyzed in the program, the results are converted into an effective temperature, which can then be used in the Cartesian farfield analysis. In this way, the pipe effects are ’hidden’ from the main routine, which simplifies the entire analysis. The third temperature array is TW (TimeStep,Pipe Length). This array holds all the water temperatures along the pipe length at each time step. The program follows a simple process all controlled by the MAIN subroutine: • Calculate coefficients: pipe convection, moisture coefficients, soil properties • Initialize arrays: water temperature, radial nodes, Cartesian nodes • Loop through all simulation days and perform the following tasks: – Update side boundary soil temperature as a function of depth, time, etc. – Update surface temperature and deep ground surface temperature. – Retrieve cycle ontime for current period in order to provide a stopping point for the simulation. – Retrieve entering pipe water temperature. – Call routine to simulate nearpipe radial system and produce an ’effective’ pipe node temperature for the farfield Cartesian system. 75 – Call routine to develop equation coefficients based on correlations and input data. – Loop through all pipe length cross sections and perform the following tasks: ∗ Loop through all nodes at the current cross section and update all nodal temperatures based on different conditions: Adiabatic on symmetric centerline, constant temperature boundaries on surface, side and deep ground This process is performed with the following required inputs: Material properties, entering water temperatures, pipe oncycle time per day, geometric data, and correlations. Some of the other key features of this simulation from a programming standpoint are the use of a function that can update a boundary temperature for both adiabatic and constant temperature boundaries. Using the same function simplifies the programming, and cleans up the code. For the entering water temperature, the same hourly boundary data was used as in the earth tube model, and is summarized in Appendix C. In Piechowski’s full program, the effects of moisture are modeled in the radial grid system. For this investigation, the effects were removed as described in section 6.2. The heart of the program is the set of algorithms used for updating the temperatures of each node at every time step. The program is designed so that at each tim 



A 

B 

C 

D 

E 

F 

I 

J 

K 

L 

O 

P 

R 

S 

T 

U 

V 

W 


