2‘. #3 {:vfs‘u‘x.) 1‘ ‘ii:-.-...- :1”, r 5 " vita ~~..~~_ ”‘7 n '.\' izgo-b gaffi‘ul #3. u}. “ . {lav 33;.iéijx. M- v ' m’i‘: u awar..§5§:v-Z : " M ' ‘.""-”‘i" u . 11"" M3. 1!"? ]‘ Afl’j » ‘4; 3:“! M. Iv" If} If}; I“, ‘ r 4.2;; "F'g '9 I” . '. .1’ . z...) 2 -..1,,'; ‘1' ”v.3f"‘/<:""/"‘«,4. f"-~' W6 '3', I.‘ ”57"”?! ’2‘?" l gr a). {erv'aflfil/f I .2131wa J .A _.! I" [f4]: ,0 .wz‘i. .‘Q‘rm’ (3".1/3 9"),ng ’ if”! W-(L ’ M a ‘ 10$ \ . #v : ‘ "1W; fiéfz’ff‘gég ‘3? .I , , Iii: _ 53.5 It ‘ k L ’ w; . «3erva ’3' mg. «.231. 1:. “3:4“ 4;» { ,5 IE: '7. “ We: JX — . ”fig; yfié’f’m Jean. aft—“4x Sflw 5%.: fig; {.15 ‘firmw ~.. . 1 .l‘h‘ifé;cfl“‘ unwensrrv LIBRARIES llllllllllllllllll‘lll l \ mm W C‘\ R 31293 007 86 4956 This is to certify that the dissertation entitled Numerical Study Of An Underground Heat Tube presented by Fauziah Sulaiman has been accepted towards fulfillment of the requirements for Ph . D . degree in Mechanical Engineering V. V Ma10r\p%fessor Date 107//é,/89 MS U is an Affirmative Action/Equal Opportunity Institution 0-12771 LIBRARY Michigan State University Ix J l 2 311/9 / 05’” PLACE II RETURN BOX to remove this checkout from your record. TO AVOID FINES rerun on or before date duo. DATE DUE DATE DUE DATE DUE _" 4:7 ——'—— :l :3 __F___: MSU Is An Affirmdivc ActhNEquul Opportunity Inflation ammo! NUMERICAL STUDY OF AN UNDERGROUND HEAT TUBE BY Fauziah Sulaiman A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Mechanical Engineering 1989 ©043blb ABSTRACT NUMERICAL STUDY OF AN UNDERGROUND HEAT TUBE BY Fauziah Sulaiman The energy, consumption of the air-to-air heat pump can be reduced, especially in winter, by using the soil as a heat source. A system of buried tube through which air is passed, has great potential in suppling higher temperature air than the ambient air, to the outside heat exchanger of the heat pump. Heat transfer from the soil to the tube, including the possibility of formation of ice lenses around the tube, was investigated over a period of time in a cold season. Models of ice formation were developed in two types of tube, the circular tube and the square tube. Latent heat released due to the formation of ice 'were included in the models. Computer simulation utilizing finite difference equations were developed, using the explicit method, where forward differences were used in time and central differences were used in space. The numerical results show the effects of increasing the moisture content of the soil, increasing the air flow rate in the tube, and the release of latent heat when soil freezes. The formation of ice around the tube played a significant role in achieving a stabilized output air temperature at a short time. Sen Scie made ACKNOWLEDGEMENTS I would like to sincerely thank Dr. Merle C. Potter, my major advisor and Chairman of the doctoral committee, for the expert guidance, suggestions and knowledge he provided me throughout the research, as it made this research both a pleasing and rewarding experience. Sincere thanks go to my doctoral committee, Dr. J. Beck, Dr. G. Ludden, and Dr. C. Somerton for their suggestions, time and support throughout this research. Special thanks go to my husband and children, for their support, patience, and prayers throughout the entirety of my doctoral study. My family in Malaysia, especially my mother, deserves special recognition, for their support and prayers. Finally, I would like to thank the Public Services Department of Malaysia and the University of Science Malaysia for providing me with a scholarship that made it possible for me to accomplish my doctoral degree. LIST LIST NClfi mm. mu. TABLE OF CONTENTS PAGE LIST OFTABLES ......OOOOOOOOOOO OOOOOOOOOOOO 0.... iv LIST OF FIGURES ......OOOOOOOOOOOOOOOOO.......... v NOMENCLATURE ... ....... ........ ..... ............. viii CHAPTER 1 INTRODUCTION .. ...................... .. 1 1.1 Background ....................... 1 1.2 Water Systems .................... 2 1.3 A Closed-Loop, Ground-Coupled system ......OOOOOOOOOOIOOOOO 3 1.4 Crawl Space Air System ........... 9 1.5 Buried Tube Air System ......... .. 11 2 THEORETICAL DEVELOPMENT .. ........ ..... 14 2.1 Introduction .......... ........... 14 2.2 Governing Heat Transfer Equations ......OOOOOOOOOOOOO 14 2.3 Finite Difference Methods ........ 18 2.4 Boundary Conditions ..... ..... .... 19 2.5 Convective Heat Transfer . coeffiCient ......OOOOOOOOOOO 22 2.6 Thermal Properties of the Soil ... 25 2.7 Initial Soil Temperature Distribution .. .......... .... 26 3 NUMERICAL MODEL ....................... 28 3.1 Introduction .............. ....... 28 3.2 Basic Equations for the One-Dimensional Model ....... 28 3.3 Basic Equations for the Two-Dimensional Model ..... .. 37 3.4 Numerical Stability .............. 62 3.5 Determination of Material Properties ......... ......... 65 3.6 Program Layout ...... ............. 66 3.7 Temperature Adjustments ...... .... 67 4 RESULTS AND DISCUSSION ....... 69 4.1 BaCkground ......OOOOIOOOOCOOOOOOO 69 4.2 Results for the One-Dimensional "Odel ......OOOOOOOOOOOOOOOOO 69 4.3 Results for the Two-Dimensional MOdel ......OOOOIOOOOOOOOOOOO 80 4.4 Comparisons for the One and Two-Dimensional Models ...... 86 5 CONCLUSIONS AND RECOMMENDATIONS 91 5.1 conCIuSions 0.00.0.0... ....... O... 91 5.2 Recommendations .................. 93 BIBLIOGRAPHY .................. .................. 95 APPENDICES ......OCOOOOOOOOOOOOOOOO 00000000000000 98 Appendix A: Computer Program for One-Dimensional Model ....... 98 Appendix B: Computer Program for Two-Dimensional Model ....... 110 Appendix C: Derivation of Heat Transfer Equations ................... 147 iii LIST OF TABLES TABLE PAGE 1 Comparisons of the Exit Temperature of the Circular Tube with Different Time Steps .... 63 2 Comparisons of the Exit Temperature of the Circular Tube with Different Space Steps inR-DireCtion ......OOOOOOIOOOOOOOOOOO...O. 64 3 Comparisons of the Exit Temperature of the Circular Tube with Different Space Steps in X-Direction O......OOOOOOOOOOOOOOOOI.0... 65 4 Values of the Soil Properties .............. 66 iv LIST OF FIGURES FIGURE 1 A Water Reinjection System ........ ..... .... 2 A Horizontal Closed-Loop, Ground-Coupled system 0.00.0.0.........OOOOOOOOO...... 3 A Vertical Closed-Loop, Ground-Coupled system ......OOOOOOOOOOOIOO... ......... 4 A Crawl Space Air System ...... ............. 5 A Proposed Buried Tube Air System .......... 6 Top View of the Proposed Two-Dimensional Tube ......OOOOOO0.00.00.00.00.0.00.... 7 Side View of the Proposed Two-Dimensional Tube 00.0.0000.........OOOOOOOOOO0.0... 8 Plane of Symmetry for the Two-Dimensional Tube IOOOOOCCOOOO......OOOOOOO ......... 9 Insulated Condition at Both Ends Of theme .....OOOOOOOOOOOOOOOO...... 10 Soil Temperature vs Depth for Michigan 8011 0.00.00.00.00.......OOOOOO....0... 11 Geometry for the One-Dimensional Model .... 12 Geometry for the Two-Dimensional Model .... 13 Node at the Center of the Circular Tube ... 14 Node Bordering the Circular Tube .......... 15 Ice Formed less than Ar/2 for the Circular we C......OOCOOOOOO0.00.00.00.00.0... 16 Ice Formed greater than Ar/2 for the Circular Tube ....-0.0.0........OOOOOOOOOOOOOOOO PAGE 10 12 15 16 23 24 27 29 30 31 33 35 36 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 Node at the Center of the Square Tube ..... Node Bordering the Center of Bottom Face of the Square Tube .................... Ice Formed less than Az/2 below the Square Tube 00.0.0000.........OOOOOOIOOOOOIOOO Ice Formed greater than Az/2 below the Square Tube 0.0.0.0000.........OOOOOOOO0...... Node Bordering the Center of Top Face of the Square Tube ....................... Ice Formed less than Az/2 above the Square Tube ......OOOOOOOOOOO0.00.0.0......... Ice Formed greater than Az/2 above the Square Tube ......OOOOOOOOOOOOOOOOOO0.0.000... Node Bordering the Bottom Corner of the Square Tube ........................... Ice Formed less than Az/2 at the Bottom Corner of the Square Tube ...... ....... Ice Formed greater than Az/2 below the Bottom Corner of the Square Tube ............. Ice Formed greater than Az/2 diagonal to the Bottom Corner of the Square Tube ...... Ice Formed greater than Az/2 adjacent to the Bottom Corner of the Square Tube ...... Node Bordering the Center of Side face of the Square Tube ....................... Ice Formed less than Az/2 adjacent to the Square Tube with Condition 1 ......... Ice Formed less than Az/2 adjacent to the Square Tube with Condition 2 ......... Ice Formed greater than Az/2 adjacent to the square Tube 000...........OOOOOOOOOOOOO Node Bordering the Ground Surface ......... Temperature of Water vs Time in Cooling Process 0.0.0.0....-......OOOOOOOOOOOOO vi 38 40 41 43 44 45 46 47 49 50 52 53 55 56 58 59 60 68 35 36 37 38 39 4O 41 42 43 44 45 46 47 48 49 50 Output Air Temperature vs Time for 10% and 50% Moisture Contents ................. Output Air Temperature vs Time for 10% and 50% Moisture Contents ................. Output Air Temperature vs Time for Air Velocities of 5 and 10 ft/s ....... Output Air Temperature vs Time for Air Velocities of 5 and 10 ft/s ...... Air Temperature in the Tube vs X-Location at 24thHour ......OOOOOOOOOOO......... Air Temperature in the Tube vs X-Location at 24thHour 0.0.0..........OOOOOOOOOOO Air Temperature in the Tube vs X-Location at 24thH°ur 0..........OOOOOOOOOOOOOOO Air Temperature in the Tube vs x-Location at 24thH°ur ..........OOOOOOOOOOOOOOOO Thickness of Ice vs X-location at 10% ”Disture Content ......OOOOOOOOOOOOOOOO Thickness of Ice vs X-Location at 50% Moisture Content ................ ...... Thickness of Ice vs X-Location at 10% Moisture Content ...................... Thickness of Ice vs X-Location at 50% ”Disture content ......OOOOCOOOOOOOO... Output Air Temperature vs Time for 10% and 50% Moisture Contents ................. Air Temperature in the Tube vs x-Location at the First Hour ......OOOOOOOOOOOOOOO Output Air Temperature vs Time for 1-D and 2-DTubes 0.0.0.0.........OOOOOOIOOOOOO Output Air Temperature vs Time for 1-D and 2-DTubes ......OOOOOOOOOOOOOOOO0...... vii 71 72 73 74 76 77 78 79 81 82 83 84 85 87 88 89 icel hicez SUI St‘W ice Ra Re NOMENCLATURE Thermal diffusivity Heat capacity Heat transfer coeeficient from soil to air in tube Ice thickness above tube in 2-D Ice thickness below tube in 2-D Heat transfer’ coefficient from. soil to air at ground surface Thermal conductivity Latent heat of fusion of water Moisture content of soil as a percentage by weight mass Nusselt number Density Prandtl number Radial coordinate Ice thickness in 1-D Radius of ice from center of tube Rayleigh number Reynolds number Temperature Time Velocity of air in the tube viii Cartesian x coordinate Cartesian y coordinate Cartesian z coordinate Change Partial differential Subscripts a atm END 2 5111' Air Atmosphere Diameter Farfield Length Soil Lo¢ation at ground surface Array Indices i x direction r direction in 1-D, y direction in 2-D time in 1-D, z direction in 2-D time in 2-D ix CHAPTER 1 INTRODUCTION 1 . 1 Background Due to the rising costs of gas and oil compared to the cost of electricity, many air-to-air heat pumps have been installed for residential use. These heat pumps are used for both cooling and heating. Since electricity is used to operate the air-to-air heat pumps, a leveling in the electric utilities throughout the year shows great promise if electricity can be used for heating as well as cooling. However, both the heat output of an air-to-air heat pump and its coefficient of performance are draStically reduced during operation in very cold weather. Thus, the use of heat pumps in regions where the heating load dominates the cooling load has been slow ,to develop. In order to improve the performance of a heat pump, various modifications were made to the heat exchanger. According to a study by the Oak Ridge National Laboratory (1), . the performance of a heat pump is significantly improved when the size of the heat exchangers is increased and the configuration of the heat exchangers is changed. These configuration changes include the use of larger evaporator flow rates, smaller compressor displacement, 1 2 smaller motor size for the compressor and the fans and finally, an increase in the number of refrigerant circuits. Installation of more efficient fans and increasing the compressor efficiency also contributes to the improvement in the performance of the heat pump. The performance of the air-to-air heat pump is primarily dependent upon outdoor air conditions. Many improvements have been made for operation in moderate climates. However, the problems caused by the operation of a heat pump in cold weather have never been solved economically, even with the improvements mentioned above. Since the air temperature is very low during the winter months, especially in northern climates, the possibility of -providing higher temperature air to the outside heat exchanger of the air-to-air heat pump is considered. This would result in higher efficiency in the performance of the heat pump and thus, would significantly increase the heating load. 1.2 Water Systems Many types of heat pumps have been tested to provide a higher temperature thermal reservoir. Surface water and well water systems have been used with water-to-air heat pumps to provide home heating and cooling (2). In winter, ground water is warmer, than air and in summer, it is cooler, Therefore, ground water acts as a heat source in the winter and a heat sink in the summer. However, a well must be drilled to tap the water so that it can be provided to the heat exchanger. A system installed by a physics professor at Ohio State University for his home, has been in use for over twenty-two years without many problems (3). It was initially meant as an experiment but after achieving the system's coefficient of performance of 3.5, it 'became a money saver. Another ground water system installed at Bear Creek, Pennsylvania, has been able to supply hot water for a 2900-square-foot house and to heat a swimming pool (4). There are some associated problems with the use of the systems, however, making them economically unattractive. The quality and the availability of groundwater ‘makes the use of this system limited (5). Ground water causes scaling problems to the heat exchanger and the possibility of aquifer depletion is always threatening. Therefore, another well is probably needed to 4 return the water to the aquifer in order to overcome the problem of aquifer depletion. This system is also known as the water reinjection system (see Figure 1). This system thus becomes uneconomical due to the high drilling costs of the wells and the restrictions due to the possibility of ground water contamination when reinjection is used. 1.3 Closed-Loop, Ground-Coupled Systems Closed-loop ground-coupled systems are divided into two basic types, horizontal and vertical systems (see HEAI EXCHANGER WATER GROUND IEVEL REINJECT ION WAT ER SUPPLY Figure l. A Water Reinjection'System. 5 Figures 2 and 3). A study (6) completed at Oklahoma State University has shown that a high coefficient-of—performance can be obtained by using the ground as a source or sink for a water-to-air heat pump. This is due to the relatively cool ground temperature in the summer and warm ground temperature in the winter. A study (7) has indicated that the horizontal type which was buried at a depth of 1 meter below the ground surface, has been able to successfully supply heated water to a home in Britain. The system has been able to supply 6 Btu annual demand of the residence. Another the 84 X 10 ground coupled horizontal system, buried 1.3 meters underground, was tested at Brookhaven National Laboratory (8). This system has been successful in supplying the heating demand of a residence at very low ambient temperatures without supplemental heating. The overall seasonal coefficient of performance of the system was 2.2. However, the performance of the horizontal type can be improved by utilizing several pipes and burying at greater depths. The geometry of the vertical type is U-shaped, configured as one deep well or several shallow wells. A SYstem with six shallow-well, U-tubes of unequal depth, Which was installed at Knoxville, Tennessee, has been able to achieve an annual performance factor of 2.39, thus reducing the maximum heating season peak demand (9). I 7 O REFRIGERATION COILS WATER OUTLET I u kl J HEAT J PUMP J WATER INLET RESIDENCE Figure 2. A Horizontal Closed-Loop, Ground-Coupled System. \A RESIDENCE HEAT ‘ WATER PUMP ——>OUTLET 1T] . I: I A._ _____ A—v v-r— ' ‘ " ' ‘— GROUND LEVEL WATER INLET U-SHAPED . TUBE V Figure 3. A Vertical Closed-Loop, Ground-Coupled System. 8 Another system utilizing thirty-seven shallow vertical couplings of 10 meters in length, was tested in Utby, Sweden, for determination of the Aground’s usefulness for seasonal heat storage (10). The resulting coefficient of performance of this system was approximately 3.0. Oak Ridge National Laboratory has done many studies, as well as designs, on horizontal and vertical ground-coupled heat exchangers. One study on a horizontal ground-coil heat exchanger with soil freezing (11) , indicated 'that 'the ‘total energy absorbed. by' a fluid is higher when the fluid inlet temperature is lower than 25°F. The study showed that at the fluid inlet temperature of 26.3'F, the latent heat absorbed was 30.6% and when the fluid inlet temperature was lowered to 10°F, the latent heat absorbed increased to 63.9%. Therefore, the increase in the latent heat was the primary reason for the increase in the total energy absorbed. Water is circulated between the heat pump and the horizontal or vertical tubes that make up the piping system that absorbs or rejects heat. In northern climates where freezing often occurs, a need of using a mixture of water and antifreeze is required for the heat pump operation. The possibility of refrigerant leaks has to be eliminated by using high quality refrigerant tubes which are quite costly. Also, small leaks must be repaired, and such work is difficult when the ground is frozen. 9 Because of the high initial cost of both the horizontal and vertical ground-coupled systems, and the relatively high cost of maintaining both systems, attention has been focussed in this study, on the air-to-air heat pump system . 1.4 Crawl Space Air Systems Energy stored in the soil beneath a house can easily be utilized by using a crawl space assisted heat pump. In Alcoa, Tennessee, such a system (see Figure 4) has been tested and found to reduce the peak heating load and the peak cooling load (12) . Air is first preconditioned by the crawl space soil in a recirculation mode before it is supplied to the outside air-to-air heat pump, therefore, resulting in a significant increase in the efficiency in the heat pump. ‘ The study revealed that the crawl space maintained a temperature near 30°F even when the ambient temperature dropped below 0'F, thereby reducing peak winter utility demand whose highest electric demand coincided with the coldest ambient temperatures. It was also found that a reduction of 26.3% in the peak heating load purchased energy and 14.4% in the total heating season purchased energy when the crawl space assisted heat pump was in use. Unfortunately, despite the great potential of this system in Producing higher temperature reservoir, the unavailability of crawl spaces in many homes is the major drawback. 10 I IN I y W I-IEAT PARTITION E: PUMP - 4.. AIR r'“ - OUT ‘— RECIRCULATING AIR TOP VIEW DWELLING INSULATION ‘_ :Issa; \\\ \\\\\\\W\\T\\\\\I SIDE VIEW Figure 4. A Crawl Space Air System. 11 1.5 Buried Tube Air Systems The energy consumption of the air-to-air heat pump can be significantly reduced, especially in the winter, by using the soil as a heat source. The potential of the soil as a moderator of the cold winter air temperatures has already been explored widely. Therefore, it is important to fully understand the heat transfer mechanisms in the soil in order to obtain the optimum benefits. A system of buried tubes through which air is passed can easily be connected to an air-to-air heat pump, as shown in Figure 5. This system has great potential in supplying higher temperature air than the ambient air to the outside heat exchanger of the heat pump and overcoming many of the problems encountered in the other systems discussed above (13). Heat transfer from the soil to the tube, including the formation of ice lenses around the tube, is investigated over a period of time during a cold period. Models of ice formation are developed in two types of tubes, the circular .tube and the square tube; this will be described in. detail in chapter 3. Computer simulation utilizing finite difference equations are developed, using the explicit method, where forward differences are used in time and central differences are used in Space. Fortran Seventy-Seven programming language is used to translate the simulation. Factors that are taken into account include: the effects of the soil moisture content, the changes in 12 \ AIR IN RESIDENCE . HEAT l (4 50 ft H Figure 5. A Proposed Buried Tube Air System. 13 soil thermal properties according to the temperature and moisture content of the soil, and the release of latent heat when soil freezes. Parameters which may be varied are the moisture content of the soil, the velocity of air in the tube, the convective heat transfer coefficient, the length of the tube, and the maximum run time. CHAPTER 2 THEORETICAL DEVELOPMENT 2.1 Introduction Two models were developed in order to determine the heat transfer from the soil to the cold air in the buried tube. The one-dimensional model has a circular cross-sectional area and the two-dimensional model has a square cross-sectional area. The tube is buried 5 feet in the ground. It was found that during a heating season, heat transfer does not take place beyond 7.5 feet (14). Therefore, a volume of 12.5 feet depth, 15 feet in width, with 50 feet length of the pipe was chosen from which to extract the 'heat (see Figures 6 and 7). The initial temperature distribution in the soil was determined by using a second order polynomial derived from Neumann's Theory (15) for’ the two-dimensional model. In the one-dimensional model, the initial temperature distribution in the soil was assigned to be constant throughout. 2.2 Governing Heat Transfer Equations The governing equation for the one—dimensional model is the one—dimensional transient diffusion equation 14 15 ‘-——O I ft U I so ft TUBE, X LY J. I“ I5 ft rl Figure 6. Top View of the Proposed Two-Dimensional Tube. 16 SOIL SURFACE I 5ft 2. . L “ ER TUBE I25 ft l4~ l5 ft a! i Fi ' - gure 7. End View of the Proposed Two-Dimensional Tube. aT._ 1 ar 8 T -—--a [ E EE'+'_"] (2.1) where the thermal diffusivity is represented by a = - ftz/hr (2.2) 5. pc and the other parameters are as follows: temperature - °F conductivity - Btu /(ft-hr-°F) = density - lbm/ft O 'O W 8 II = specific heat - Btu/(lbm-°F) The governing equation for the two-dimensional model is the two-dimensional transient diffusion equation 2 2 6T 6 T 6 T 3E=°‘[—§+—§] (2'3) ay 62 Heat transfer in the x-direction is quite small and has been neglected. Boundary conditions and the initial condition needed in the Solution of the above diffusion equations will be Presented in a later section. A finite-difference method, which is often used in solving partial differential equations, was employed 18 using forward differences in time and central differences in space; it is known as the explicit method. Solutions were obtained by placing a finite number of points in a selected grid pattern within the soil and air in the tube. At these points, or nodes, finite-difference equations for the temperatures were obtained by direct replacement of the partial differential equations or by an energy balance on an element. The error of the method is on the order of size of space step squared and the time step. 2.3 Finite Difference Methods There are several finite difference methods that can be used to solve the heat diffusion equations. Descriptions of some of the methods follow: 1. The Implicit Method is similar to the explicit method outlined above, but instead, uses a backward difference approximation for the time step. It has the same order of error as the explicit method but is stable unconditionally. However, this method cannot be used to solve for a single unknown in terms of previously calculated values. It was also found that the method is not efficient for multi-dimensional problems (16). 2. The Crank Nicolson Method, which is known as the 1mProved implicit method, is also unconditionally stable. This method averages the forward difference and backward 19 difference approximation for the time step. Unfortunately, this method produces a system which is no longer tridiagonal and thus, requires unnecessary computations (17). 3. Alternating Direction Implicit Method (ADI) of Peaceman and Rachford was proposed to overcome the problems encountered in the explicit and the implicit methods. This method uses an implicit method in one direction and uses an explicit method in another direction. It has, however, a bias on intermediate values and therefore is not to be greatly trusted (18). . According to a preliminary study, the purely explicit method was rated third from a comparison of nine methods of solving the heat diffusion equation (19). This was based on relative accuracy, ease of programming, computation time and computer storage requirement. However, the explicit method has a stability criterion which must not exceed 0.5. Due to this, substantial computing time will be needed since there is a need to employ small time steps. This is the method selected in this study. 2.4 Boundary and Initial Conditions The boundary conditions in the one-dimensional model are as follows: 20 1. At r = r1 (see Figure 14), =k air) 5 6r (2'4) h (Ts - T where h is the convective heat transfer coefficient from the soil element to the air in the tube (the thin wall tube material was ignored) and k5 is the soil conductivity. T s and Tair are the temperature of soil and air, respectively. 2. At r = rF (farfield), T = T (2.5) where TF is the given constant farfield temperature. The initial conditions at t = 0 are (2.6) where the initial temperature of air and soil are assigned to be the farfield temperature at all locations. Boundary conditions for the two-dimensional model are: 21 h’ h’ . 1. At y = h’ and -§ < z < 5 (see Figure 29), aT kB 5? h (Ts Tair) (2.7) h’ h’ . 2. At z = -h’ and -E < y < 2 (see Figure 18), 6T _ _ -ks 52 _ h (Ts Tair) (2'8) h’ h’ . 3. At z = h’ and -E < y < 5 (see Figure 21), aT _ _ ks 52 _ h (Ts Tair) (2'9) 4. At y = 0 and z > h’ and z < -h’, _= O (2'10) The initial conditions are: Tair = Ts(z = tube height) (2.11) T8 = Ts(z) The air inlet condition is: Tair(x = 0,t) = O F (2.12) 22 As above, h is the convective heat transfer coefficient and k5 is the soil conductivity. The atmospheric air temperature is always maintained at 0°F and the soil temperature at 12.5 feet depth is always set at the deep soil temperature. Due to symmetry, heat transfer need only be calculated on one side with a vertical plane bisecting the tube longitudinally (see Figure 8). An insulated condition is being imposed at each end of the tube since heat transfer there is assumed small and thus ignored (Figure 9). This condition is considered conservative as there is actually heat transfer occuring at each end of the tube. 2.5 Convective Heat Transfer Coefficient The velocity of air flow in the tube is assumed to be sufficiently high to ensure a turbulent flow. The turbulent flow in the circular smooth tube is considered to be fully developed. The Nusselt number is thus calculated by using the Dittus-Boelter equation (20), NuD = 0.023 Reg/5 Pr°°4 (2-13) Where ReD is the Reynolds number (based on diameter of the tube and average air velocity) and Pr is the Prandtl number. The convective heat transfer coefficient is thus obtained. F°r the square tube, the same equation is used but the 23 )— PLANE OF SYMMETRY Fi gure 8. Plane of Symmetry for the Two-Dimensional Tube. 24 Figure 9. 'Insulated Condition at Both Ends of the Tube. 25 diameter of the circular tube is replaced by the hydraulic diameter of the square tube. To be conservative, at the soil surface, free convection (worst case condition) is assumed. The correlation used for the calculation of Nusselt number is obtained by the improved correlations of McAdams (21), Nu = 0.15 Rel/3 L L (2.14) where RaL is the Rayleigh number. 2.6 Thermal Properties of the Soil The soil conductivity and thermal diffusivity are relatively low in dry soil. However, moisture content in the soil increases the thermal conductivity and thus provides for increased heat transfer (22) . In addition to moisture content, soil type and dry density must be known in order to estimate the thermal conductivity of the soil. Since there is a large variation even in similar soils, soil thermal properties were empirically determined. During the process of freezing, heat capacity of the soil is significantly high. This is due to the release 0f latent heat of fusion in the moist soil. Frozen soil has the highest thermal conductivity which is 0.833 Btu/(hr-ft-'l'-') but a relatively low heat capacity of 0.2 Btu/(1bm--'F) . The high thermal conductivity results from the soil moisture completely changed to ice and that ice is 26 more conductive than water. waever, heat capacity of ice is lower than water, thus resulting in low heat capacity. Thawed soil has relatively low thermal conductivity and heat capacity which are 0.458 Btu/(hr-ft-°F) and 0.3 Btu/(lbm-'F), respectively (23). 2.7 Initial Soil Temperature Distribution Comparison with actual data of’ Michigan soil temperatures (24) and Neumann’s theory (25), led to a linear distribution merged with a second-order polynomial for the two-dimensional model (Figure 10). Irt was noted that when the ground surface is frozen, the temperature variation in the frozen surface layer of the soil is linear, whereas in the thawed condition, the temperature variation. in. the surface layer is sinusoidal (26). However, this sinusoidal variation does not extend to a significant depth so it was not taken into account. In the one-dimensional model, the - temperature in the soil was assigned a constant value equal to the constant farfield temperature. DEPTHIFEET) 30 l2 ~ Figure 10. 27 TEMPERATURE (°F) 32 4O NEUMANN Soil Temperature vs Depth for Michigan Soil. CHAPTER 3 NUMERICAL MODEL 3.1 Introduction Two types of models were developed, one with cylindrical coordinates (the one-dimensional model) and the other with Cartesian coordinates (the two-dimensional model). The geometry for the numerical models are shown in Figures 11 and 12. The possibility of ice formation around the tube was taken into consideration and modelled appropriately. The two models provide for a comparison with the two-dimensional model being the more accurate. It is anticipated that the one-dimensional model may' be sufficiently accurate. 3.2 Basic Equations for the One-Dimensional Model The equations which describe the heat transfer in the one-dimensional model are as follows: 1. When. the node is at the center of the tube (Figure 13), the following equation derived from a steady-state energy balance is used (see Appendix C): 28 29 4r 'T—1 -=r====1 Ar Ax IADUMBATKZ SURFACES Figure 11. Geometry for the One-Dimensional Model. g7 >Y. //, AY A z I__—T 152 fiAx ADIABATIC SURFACES Figure 12. Geometry for the Two-Dimensional Model. 31 a, Figure 13. Node at the Center of the Circular Tube. 32 k _ 2hAx k _ k k Ti+1,0 _ rm ( Ti,1 Ti,0 ) + Ti,o (3°11) a pa 1 where pa is the density of air, cpa is the heat capacity of air, u is the velocity of air and r1 is the radius of the tube. The superscript "k" denotes the time and the subscripts denote x and r, respectively. 2. When the node borders on the tube (Figure 14) , the following equation is obtained with application of the boundary condition given by equation 2.4 (see Appendix C): k+1 At k k T. = h 4nr Ax (T. — T. ) 1,1 A; [ 1 1,0 1,1 psn(r1 + 4 )ArAxcps Ar + k 4n(r1 + §—)Ax (T.k - T.k ) (3 2) s2 Ar 1,2 1,1 ° k + Ti,1 where ps is the density of soil, cps is the heat capacity of soil and k82 is the thermal conductivity of soil evaluated k k . at (Ti,j+1 + Ti,j)/2 . The calculated temperature will always be checked to determine if some of the soil is k+1 < 32°F and r. < Ar/2 , part of the soil fro . zen If Ti,1 ice 33 b-X Figure 14. Node Bordering the Circular Tube. 34 element next to the tube will be frozen (Figure 15). The equation used for the ice formation is (see Appendix C), Ar. _ Ar ice - pSLM 2nR Ax At n(r1 + z—)ArAx(T At [ fgfgg k+1 _ k ] where L is the latent heat of fusion of water, M is the moisture content of the soil as a percentage by weight and R is the radius of ice formed from the center of the tube. If the element is completely frozen, the element directly next to it may freeze partially (Figure 16). The following equation is then used for the ice formation (see Appendix C): At pscps k+1 k ce psLM zfifi Ax[ 2nr ArAx (Ti, T. ] Ari At 2 2 where ksl is the thermal conductivity of soil evaluated at (T§Ij_1 + Tilj)/2 . 35 e Circular Tube. Figure 15. Ice Formed less than Ar/2 for th ._—___ j 36 Ina TFAJKG Figure 16. Ice Formed greater than Ar for the Circular Tube. 37 3. For the remaining soil elements, the following equation is used (see Appendix C): Ar T¥*? = At [ k 2n(fi §_)Ax (Tk - Tk ) 1,3 mscps sl Ar 1,3-1 1,3 2n(rj + %£)Ax k k + ksz Ar (Ti,j+1 ' Ti,j) ] (3'5) k * Tm 3.3 Basic Equations for the Two-Dimensional Model There are several equations that describe the heat transfer in the two-dimensional model. They are as follows: 1. As in the one-dimensional model, steady state is also assumed with an energy balance on the air element in the tube. The following equation is derived for the node at the center of the tube (Figure 17): l Ax [ l l T. = —————— 2h AXA . - . 1+l,0,0 macpau z (T1,-l,0 T1,0,0) + 1 1 2“ AXAZ (T1,1,0 ' T1,0 0) + 2h A A T1 1 3 5 x Y ( i,0.-1 ' Ti,0,0) ( ' ) 38 +N Figure 1?. Node at the Center of the Square Tube. 39 l l + 2h AxAy (Ti,0,1 - Ti,0,0) ] 1,0,0 where ma is the mass of air in an element; the superscript again refers to time and the three subscripts to x, y, 2, respectively. 2. When the node borders on the center of the bottom face of the tube (Figure 18), the following equation is used (see Appendix C): TiTO,-1 z EAE“[ ks1A§%E (Ti,-1,-1 ' Ti,0,-1) 5 ps + ks2AI$£ (Ti,1,-1 ' Ti,0,-1) + 2ks3é%%x (Ti,0,-2 ' Ti,0,-1) (3'7) + 2h AXAY (Ti,0,0 ' Ti,0,-1) ] + T1,0,4 1+1 If 'r. ° . 1,0,-1 < 32 F and hicez < Az/2 , the $011 element directly next to the node is partially frozen (Figure 19). The following equation is used for the ice formation (see Appendix C): 1‘; A \\\\\\ . \‘ i , \\\ “A Face of the Pi ure 18 Node Bordering the Center of Bottom 9 Square Tube. —» 4_________ j 41 FN "Y .\\\ \‘h‘. hiced Ahicez —cI —q t T—)>X Figure 19. Ice Formed less than Az/2 below the Square Tube. 42 mc Ah. = At _§_E§ (T%+1 - T; ) 1ce2 2pSLM AxAy At 1,0,-1 1,0,-1 (3.8) If the element is completely frozen, the element directly below the frozen element may be partially frOzen (Figure 20). The equation for the ice formation is as follows: m c Ah. = At -§—E§ (Ti+1 - T} ) 1ce2 pSLM AxAy At 1,0,-2 1,0,-2 (3.9) 3. When the node borders on the center of the top face of the tube (Figure 21), the equation used for heat transfer is similar to equation (3.7). The only difference is that a different boundary equation is applied whereby the soil-to-air convection term is positioned appropriately in the equation. The formation of ice is also similar but labelled differently (Figures 22 and 23). 4. When the node is at the bottom corner of the tube, special consideration is taken to account for the corner element with surface convection (Figure 24) . For no ice formation the following'equation is used (see Appendix C): 43 h- 1ce2 -- -- i AI‘icez I 2t F X h- I IceZ ‘-.------ im‘icez I Fi(lure 20. Ice Formed greater than Az/2 below the Square Tube. 44 ’N y-Y ANAL s-b-X afa— FiQUre 21. Node Bordering the Center of Top Face Square Tube. of the 45 ’N I.— l.— hIceI AhiceI P'Y 59/ .\\ AhiceI him i .I—u- ~fi§tX ->IIo— Figure 22. Ice Formed less than Az/2 above the Square Tube. 46 Z A "L- ‘i-Ahiclu hkeI I / --------- iflAhmeI “MEI T >X Figure 23. Ice Formed greater than Az/2 above Tube. the Square 47 ’N \h\\ V Figure 24. Node Bordering the Bottom Corner of the square Tube. 48 Titi,-1 = 3%AE“ [ kslé%%£ (Ti,0,-1 ’ Ti,1,-1) s ps + 2kszé%%é (Ti,2,-1 ‘ Ti,1,-1) + 2k539%%1 (Ti,1,-2 ' Ti,1,-1) + ks4é%%x (Ti,1,0 ' i,1,-1) (3'10) + h AxAz (Ti — T T 1 ) 1,0,0 i,1,-1 l i 1 + h AxAy (Ti,0,0 T ,1’_1) ] 1 + Ti,1,-1 When there is s0me ice formation occuring in the element, h < Az/2 (Figure 25), the following equation is then used ice2 to account for the release of latent heat of fusion (see Appendix C): 1+1 2At [ T. _ = -————— p LM Ax (Ay + Az 1,1, l 3mscps s Ah. 1ce2 + 4hice2 + 2Ahice2) At ] + T; (3 11) 1,1,-1 ° If . Az/z 5‘ h1ce2 frozen corner element will be partially frozen (Figure 25)- : 3Az/2 , the element directly below the The following equation for that element is: 49 —’N y-Y N‘Q Ihkez Fnggti t th :re 25 . Ice Formed l ess than A f S . z/2 a e Bottom Corner 50 pN "‘Y “' Ahic 2 F‘ - J~9ure 26. Ice Formed greater than Az/2 below the Bottom Corner of the Square Tube. 51 Ah. 1+1 At [ 1ce2 T. _ = p LM AxAy —————— J 1,1, 2 mscps s At + T1 (3.12) i,1,-2 The equation used for the element diagonal to the frozen corner element (Figure 27) is as follows: 1+1 _ At _ 1 Ti,2,-2 - fi_E__ [ psi"M Ax [ 2(hice2 2A2) 5 PS . + Ah ] Ahice2 ice2 A5 + T1.“ (3 13) 1,2,-2 ° The equation used for the element directly to the right of the frozen corner element (Figure 28) is: Ah 1+1 At [ ice2 T. _ = ————— p LM AxAz —————- ] 1,2, 1 mscps 3 At + T1 (3.14) i,2,-1 52 Phi A \N 4 Y 7 Ahicez Jpn.- NIcez Figure 27. Ice Formed greater than Az/2 diagonal to the Bottom Corner of the Square Tube. 53 'N \\Q R 7 —-Il‘— ' Y I AI‘icez I? hice2 Figure 28. Ice Formed greater than Az/2 adjacent to the Bottom Corner of the Square Tube. 54 5., When the node is at the top corner of the tube, similar equations are used for the frozen and unfrozen states. For the nodes surrounding the frozen top corner element, similar equations are used for the heat transfer in the element, taking into account the release of latent heat of fusion. 6. When the node borders the center of the side face of the tube (Figure 29) , the following equation is used when no ice formation occurs: 1+1 _ At 1 _ 1 T1,1,0 ‘ mscps [ 2h AXAZ (Ti,o.o Ti,1.0) AxAz l l + stz Ay” (Ti,2,0 ' Ti,1,o) AxAy 1 _ 1 + ks3 Az (Ti,1,-1 Ti,1,o) (3°15) AxA l l + ks4_32[ (Ti,1,1 ‘ Ti,1,o) ] l + T1,1,0 v01: . ' en 0 < hicel s Az/2 and hicez 0 (see F1gure 30), the following equation is used: 1+1 a Ah. + Ah. '1‘. = At 2 2 1ce1 1ce2 1.1.0 mscps [ Pam A" [(A2’ + hicel At ] 1 + . . Th1,o (3 16) 55 —pN I. T PY Figure 29. Node Bordering the Center of Side Face of the Square Tube. 56 ‘T‘iIfl hiceI _ r Y // (N‘iceI + Ahicez) / 7- \\\ IPUigndre 3 0. Ice Formed less than Az/ ' . 2 ad acent to t Tube w1th Condition 1. 3 he Square 57 When 0 < h. s Az/2 and 0 < h. 1cel 1ce2 Az/2 (see Figure 31),the following equation is then used: 1+1 _ At / 2 _ 2T Ti,1,0 — mscps [ psLM Ax (A2) + (hicel hicez) Ahicel + Ahice2 At + T]: (3 17) 1,1,0 ° When the element is completely frozen, the next element directly to the right of the frozen element may be partially frozen (Figure 32). The equation used for the heat transfer in that element is: 1+1 At 1 2 2T T. = __ - .. 1.2.0 mscps [ 2951"M Ax ¢r(Az) + (hicel hicez) Ahicel + Ahice2 At + T1 (3 18) i,2,0 ° 7. When the node borders the ground surface (Figure 33). the following heat transfer equation is used: 58 / "— hicel /r I / {.— hice2 A’Y \N Ice Formed less than Az/2 adjacent to the Square Figure 31 . Tube with Condition 2. 59 ’N hiceI L\\\ “ken Figure 32. Ice Formed greater than Az/2 adjacent to the Square Tube. Ly§;$)\\‘ \ 1\ \\\\ R“ \\\N y >Y //, J////// // \ axe— \ \ \ 2L 1_’,x Figure 33. Node Bordering the Ground Surface. 61 Tifzjasur = 5%; I: kslAjfi—z mi,j-1,zsur - Ti,j,zsur) + kszA—ng (Ti,j+1,zsur - Ti,j,zsur) + 2ks3A—)A‘_gx (Ti,j,zsur-1 - Ti,j,zsur) + 2hsur AxAy (Tatm - II.i,j,zsur) ] + i,j,zsur (3'19) '1 where hsur is the convective heat transfer coefficient of the soil to the outside air. 8. For the other soil elements, the following equation is used (see Appendix C): Tifj,k ’ EA§" [ £121?E (Ti,j-1,k ' Ti,j,k) sps + kszéi$5 (Ti,j+1,k ’ Ti,j,k) + “539%?! (Ti,j,k-1 ‘ Ti,j,k) (3°20) + ks4é%%x (Ti,j,k+1 ’ Ti,j,k) ] + T1 j'k Where ksl is the conductivity of soil evaluated at J. (Ti,j-1,k + Ti,j,k)/2 , k$2 is the conductivity of soil 62 1 i,j,k)/2 ' ks3 /2 and finally, RS4 is . . . l l the conduct1v1ty of $011 evaluated at (Ti,j,k+l + Ti,j,k)/2° l . . . evaluated at (Ti,j+l,k + T is the conduct1v1ty . l l of $011 evaluated (Ti,j,k-1 + Ti,j,k) 3.4 Numerical Stability The numerical model which was developed using the explicit method has a stability criterion. For a stable numerical solution, the formula used for the one-dimensional model is —EAE§ s 1% (3.21) (Ar) and for the two-dimensional model, the formula used is «At 5 (A102 + (Az)2 (3.22) Ana where Ay = A2. For reasonable accuracy, the time step and the space steps have to be made quite small (27). Several tests were conducted to select the appropriate time step and spaCe step for the one-dimensional mOdel considering only the exit temperature of the tube. Tables 1, 2 and 3 indicate the tests that was conducted 311d selections that was made. It was found that a time step of 10 minutes and a space step of 5 ft in the x-direction would give reasonably accurate results. The space step in the 63 r-direction is varied between 0.2 ft and 0.4 ft, allowing only 2 soil elements to freeze. Table 1. Comparisons of the Exit Temperature of the Circular Tube with Different Time Steps. Time At = 1 min At = 10 mins At = 1 hr Time Tout Tout out (hr) ('F) ('F) (°F) 0.5 14.7540302 15.2320642 1.0 13.5392723 13.7052917 16.5444374 1.5 13.4371080 13.4371080 2.0 13.4371080 13.4371080 13.4371080 2.5 13.4371080 13.4371080 3.0 13.4371080 13.4371080 13.4371080 3.5 13.4371080 13.4371080 4.0 13.4371080 13.4371080 13.4371080 4.5 13.4371080 13.4371080 5.0 13.4301677 13.2995367 12.4048958 5.5 12.9083624 12.8521147 6.0 12.4164248 12.3375645 11.0715809 Here, Ar = 0.564 ft and Ax = 1.0 ft. The selected value for At is 10 mins. 64 Table 2. Comparisons of the Exit Temperature of the Circular Tube with Different Space Steps in the x-direction. Ax = 1 ft Ax = 5 ft Ax = 10 ft Time Irout Tout Tout (hr) ('F) (°F) (’F) 0.5 15.2320642 15.4745922 15.8068676 1.0 13.7052917 13.9133835 14.1865015 1.5 13.4371080 13.6633844 13.9620523 2.0 13.4371080 13.6633844 13.9620523 2.5 13.4371080 13.6633844 13.9620523 3.0 13.4371080 13.6633844 13.9620523 3.5 13.4371080 13.6633844 13.9620532 4.0 13.4371080 13.6633844 13.9620532 4.5 13.4371080 13.6633844 13.9620532 5.0 13.2995367 13.5441265 13.9620532 5.5 12.8521147 13.0398951 13.4373751 6.0 12.3375645 12.4092398 12.8577614 Here, At = 10 mins and Ar = 0.564 ft. The selected value for Ax is 5 ft. Table 3. Comparisons Ar = 0.2 ft 65 Ar = 0.4 ft of the Exit Temperature of the Circular Tube with Different Space Steps in the r-direction. Ar =0.564 ft Time Tout Tout Tout (hr) ('F) ('F) (‘F) 0.5 13.6781063 14.8539438 15.4745922 1.0 13.6633844 13.6633844 13.9133835 1.5 13.6633844 13.6633844 13.6633844 2.0 12.9093580 13.6633844 13.6633844 2.5 12.0458002 13.6633844 13.6633844 3.0 11.3234854 13.6633844 13.6633844 3.5 10.8421450 13.4894409 13.6633844 4.0 10.5065155 12.7473459 13.6633844 4.5 10.3096199 12.0335646 13.6633844 5.0 10.2087717 11.3665028 13.5441265 5.5 10.1735430 10.6625042 13.0398951 6.0 10.1728191 10.0137520 12.4092398 Here, At = 10 mins and Ax = 5 ft. The selected value for Ar is between 0.2 and 0.4 ft, allowing only 2 soil elements to freeze. 3.5 Material Properties Soil properties 'vary enormously and. therefore empirical values must be determined (28) . These values are tabulated in Table 4 ' 66 Table 4. Values of the Soil Properties. Normal Soil Conductivity 0.458 Btu/(hr-ft-°F) Frozen Soil Conductivity 0.833 Btu/(hr-ft-‘F) Normal Soil Heat Capacity 0.3 Btu/(lbm-°F) Frozen Soil Heat Capacity 0.2 Btu/(lbm-°F) Soil Density 95.0 lbm/ft3 The air outside which is being fed into the tube is selected as a winter condition of 0'F at all times. Since the air in the tube may vary between O'F to 40'F, the density of air at 20'F is chosen which is 0.08275 lbm/ft3(29). Atmospheric pressure is used since pressure varies only slightly from atmosphere. The corresponding heat capacity of air is 0.24 Btu/(lbm-'F). The latent heat of fusion of water is 143.3 Btu/(1bm-'F) and the heat capacity of water is 1.0 Btu/(lbm-'F). The heat capacity of ice is 0.5 Btu/(lbm-°F). 3.6 Program Layout Initially, some parameters are input interactively so that several runs can. be made without changing the contents of the program. The parameters are the element size, the moisture content in the soil, the velocity of air in the tube, the stability constraint and 67 the convective heat transfer coefficient from the soil to the tube. The initial temperature distribution in the soil is then calculated and the temperature of air in the tube is assigned to be the temperature of soil at that depth. The program next selects the heat capacity and the thermal conductivity of soil according to the temperature surrounding the node. Then it will select the appropriate heat transfer coefficient for calculation of the temperature at that node. Calculation of temperature of the air in the tube is accomplished only after the soil temperature is found at that time step. 3.7 Temperature Adjustments Ice formation. around the tube occurs if the temperature of the nodes bordering the tube is less than 32’F. Temperature is always reset to 32’1" if the soil element is not completely frozen. When the element is completely frozen, the temperature will be calculated and will be less than 32'F. According to Jumikis (30), latent heat of fusion is released as long as there is water left in the soil element in the process of cooling (Figure 34). The release of latent. heat. of fusion. ‘will contribute Significantly to the energy transfered to heat the air in the tube (31). TEMPERATURE 68 TIME WATER —_—-| WATER AND ICE I-——— SOLID ICE Cooling ——-> Release of latent heat of fusion FiGure 34. Temperature of Water vs Time in Cooling Process. CHAPTER 4 RESULTS AND DISCUSSION 4.1 Background Several parameters were varied in order to obtain data for comparisons. The parameters that were varied are the air flow rates, the convective heat transfer coefficient at the surface of the tube, and the moisture content in the soil. For the one-dimensional model, the computer simulation was in operation for 24 hours. However, considerable computing time is needed to run the two-dimensional model; and therefore, a run only between 1 to 2 hours was done to check the exit temperature with the one-dimensional model. Since the models were developed to simulate operation during a cold winter period, the air inlet temperature is set at 0°F. Criterion for termination of the program is when the maximum time of continuous operation is done. For all the one-dimensional model runs, a maximum time of continuous operation of 24 hours was set. 4-2 Results for the One-Dimensional Model The simulation predicted that a temperature of 12-4’F is obtained for 10 percent moisture content in the 69 70 soil and 13.7'F for 50 percent moisture content in the soil at an air flow rate of 300 cubic feet per minute (u = 5 ft/s) after 24 hours of continous operation as indicated in Figure 35. This figure is a graph of the output air temperature of a 50 ft circular, 6.77-in-radius tube at an air velocity of 5 ft/s for different moisture contents in the soil as a function of time. Here, the output air temperature difference after 24 hours is 1.3°F. Figure 36 shows the output air temperature of the 50 ft tube as a function of time for different moisture contents of the soil at an air flow rate of 600 cubic feet per minute (u = 10 ft/s). At 10 percent moisture content of the soil, the output air temperature of the tube is 8.5°F and at 50 percent moisture content of the soil, the output air temperature of the tube is 10.2°F after 24 hours of continuous operation. In this case, the difference in the output air temperature is 1.7’F. For both cases of 5 ft/s and 10 ft/s air velocity, the output air temperature of the tube for 50 percent moisture content in the soil is higher than the output air temperature for the 10 percent moisture. This is due to the higher heat capacity of the soil which depends on the soil moisture. Figures 37 and 38 are the graphs of output air temperature of air velocities of 5 ft/s and 10 ft/s as a function of time at 10 and 50 percent moisture contents in the soil, respectively; At the beginning, the output temperature difference between the 5 ft/s air and 10 ft/S 71 +N ( L . .L x l \ A. _ 1‘ Q _ 3505 9.5 mu 3 m. o. .2 N. o. a o .v . . _ . . . q . . . woa won N. n. .v. m. m\: m m. on: c. > p :oo_m> :< can; £95.. c. : on lo we E 85686.0-ch o— )GJHLEJGCLUGL 1W 1nd1no :! ( 72 3.505 we; +N Nu on m. 2 t N. o. a o _ . . . . . . . . _ «0H mom 0. N. m. +_ 1W 1nd1fl0 8 ) GJnLBJedLU .J ...: 73 +N .2: 2 tan. a .3 3 .‘ " l '4!— :oo_m> =< .2 we; m> @5659: me =< 39:0 Km m 59“. NN 6:55 we: 9.. m. w. +_ .2 o. m o e N . . . mxum oH mxuu m . . o— N. m. ¢‘ m. w. (3.)em121edwei NV 1nd1nQ 74 .2: 2 was m .o mm;_oo_v> :< lo. 9:: 9 053.882 :< 59:0 .mm 0.59”. Amezoé 9:: ...N «N on m. 2 .I N. o. m o . .v . . . . . . . . . . . m\uu ca m\uu m O. N. n. +_ m. 2 .\7' NWde edwei Jl (:4.) 61mm 75 air was approximately 1.5'F for both cases. After about 12 hours of continuous operation, the difference became close to 3.5°F for the 50 percent moisture content in the soil (see Figure 38), and 3.8°F for the 10 percent moisture content in the soil (see Figure 37) . For the one-dimensional tube of 50 ft in length with air velocities in the tube of 5 ft/s and 10 ft/s, Figures 39 and 40 show graphs of air temperature in the tube as a function of x-location along the tube at the 24th hour of continuous operation. It appears that for different moisture contents of the soil, the air temperature in the tube does not vary significantly at each x-location along the tube. For both cases of the air velocity in the tube of 5 ft/s (see Figure 39) and 10 ft/s (see Figure 40), the temperature difference between the 10 and 50 percent moisture content, does not exceed more than 1.7°F throughout the length of the tube. Air temperature in the tube versus x-location along the tube at air velocities of 5 ft/s and 10 ft/s for 10 and 50 percent moisture contents in the soil is depicted in Figures 41 and 42, respectively. For the 10 percent moisture content in the soil, the temperature difference for the two air velocities along the tube varies between 0.54“}? and 3.8'F and for the 50 percent moisture content in the 5011, the temperature difference varies between 0. 39'F and 3.51'F. 76 .50: £2 E cozqoobk m> .33 on. E 228883 __< .mm 039“. 2%: :o_§o._-x ”0.." “on m3 m m. 23 c. 3.oo_m> :< con; 59.3 c. : Om so we: .mco_mcme_oumco O— 2 (i) 63an at“ U! GJMBJSdUJ-ei ng 77 .So: Evin co_.moo-Tx m> we; 9: c. moiflmaefl =< .0: 2:9“. cmmc COZMQOJIX MN.» .93. mfim m..~m MNN m.- .92 ”we . ms mN . . _ _ q . . . fl . «on i ll m3 9 m. we: :. >:oo_o> s< can; 59.3 c. : 0m .0 one: 6.5.2065-ch 0. 5 (i) ‘3an 3‘41 U! eJnLEJedwei ng 78 «Kw m.N¢ .50: Svfim cozmoo-Tx m> one; 9.: E 222863 :< J? Ezoi 2mm: cosmoomnx mxm MN... 93 m.- m5 mi 4 ms .m..N M\um o." m\uu m . . . . . _ . 22:00 22922 $0. :m £95.. 5 : on .o .33 6:06:35:ch O. (2L) 53an 3H1 U! aanJadwei ng 79 .So: Evin co._moo._|x m> we: of E 222888 :< .Nv 2:9“. :8: 8:83; Wk. 9% hem m.~m msm 9.3 m5 .....N. _ ms ma _ . . . s . . . 4 '1': 0 m .. m who 3 . . 1 o. 4 who. m if Iii Ill m. 22:00 @2202 $8 6 595.. c. : on .o 83 585885-30 (i) eqni GUI U! eJmeJe-dwei my 80 Figures 43, 44, 45 and 46 show the thickness of ice along the length of the tube at different air velocities and moisture contents of the soil. At an air velocity of 5 ft/s, Figures 43 and 44 show the ice formation at the 8th, 16th and 24th hour at 10 and 50 percent moisture content in the soil, respectively. Figures 45 and 46 show the thickness of ice formation at the 8th, 16th and 24th hour at an air velocity of 10 ft/s with 10 and 50 percent moisture content in the soil, respectively. It is observed that there was more ice formed when the moisture content of the soil is 10 percent for both cases of the air velocity. This is due to more time needed to freeze the moisture in the soil for higher moisture content. 4.3 Results for the Two-Dimensional Model Criterion for termination of the two-dimensional model program is when the output air temperature achieved do not change from one time step to the next. This stability '10 0F. criterion is set at 1.0 x 10 Figure 47 shows the output air temperature as a function of time for 10 and 50 percent moisture content of the soil when the air flow rate is 600 cubic feet per minute. For the 10 percent moisture content in the soil, a stabilized output air temperature of 14.96‘1? was obtained after 1.09 hours of continuous operation whereas for the 50 Percent moisture content in the soil, a stabilized output air temperature of 14.98'F was obtained after 1.98 hours of 81 .3250 $3902 $9 .6 co_.aoo.Tx m> m2 .0 3053:: .mv 239m 2%: cosmoofx Wt. m..~¢ mém ”NM. mKN m..NN m5 muN. . ms. mm _ u . . m . . . . . . O u 1 t... m 1 ..o .... ... i 1 2 Wm. N6 1 m6 1 i P .2 4: m a 2.: 5 rooms =< 8;; 59.3 s : cm .0 we: .mcoacmsélmco ( 1925) 331 ;o ss'euwolui 82 ._._N.___:.,_.. 3:..........s. 30m :. ......Bc._-x a, 8. .o 395:: .3 9.5; :8: 8:83; 95¢ .98 9%. m.~m msm mdu m5 m8 . m5 m.~ .1]. . . . . _ _ J _ 7 . LJ Q— 7 1 t. +N h L i 1h"! m). m m. 63:-.. :_ >__.,5_o> _.< :9? 595-. c. : Om .o 2:: .mcoacosfioco ..0 N6 83 ....ona. 23222 80. 2. ..o:so3-.x .9 oo_ .o mmmcxoip .mg 239“. 28: 8:80-; m? mi» MR mum TR n.3, m5 m8: . .95 m8 0 a . . _ . . . q _ . g ..~.o E’s . 1].! A H... fr / 1 .... +~ 1.8 i 8.0 m2. 0. m. 33 E >:oo_o> .2 cos; 595.. c. : On .0 we: .mco_mco§o-oco 84 maniv 4:2:ch 23922 wom E :o._831x ms oo. .o 3936:: .ov 939.1 cwmb COZMUOITX h.~+ mfim. 9N». mfiw .méN .95 h.~. , m6 m.N . . . . . . _ . _ .5 2 x: +.N 1 so m2. 9 m. 2.3. c. >:oo_o> .2 cos; 59.3 E : 0m .o 33 .mcoficoczonocO No 4 85 ..cmEau 0.5.90.2 $5. .5... ago. 5. 05; m» 3:.309:o.. .7... 50.30 He 05?“. $505. 0E; ww mam Wm 04 we m.P 0a 50 v.0 ..o _ . . . n . . . . . . . . . . . . . . 05.20:. 05.905 cson... .x. or mi 0.9 09 0.0.. m... o. m. 83 :. >..00_m> ..< :22, £95.. 5 .. 0m .0 we: 38.308502; (:1.) .9an19de m1 1nd1n0 contir simulz tempe: percen diffe: the t in th« depiC' diffe: x-loc. conte: varia quick obser 4.4 as a are g CaSes minut. perCel is 1 tw0‘dj 86 continuous operation. After the beginning of the simulation, at 0.5 hours of the tube's operation, the output temperature difference between the two cases was 4.1 percent. However, after stability, the temperature difference was only 0.2 percent. Air temperature in the tube along the length of the tube is compared for 10 and 50 percent moisture content in the soil when the air velocity in the tube is 10 ft/s, as depicted in Figure 48. As in the one-dimensional tube, the difference in the air temperature in the tube at each x-location along the tube for the different moisture contents in the soil, does not vary significantly. The variation is only between 0.04°F and O.2°F. Since the output air temperature stabilized very quickly for both runs, the ice formations could not be observed due to program termination. 4.4 Comparisons for the One and Two-Dimensional Models Output air temperature in the 50 feet long tube as a function of time for the one and two-dimensional models are graphed in Figures 49 and 50 for comparisons. For both cases, the air flow rate in the tube is 600 cubic feet per minute. Two moisture contents of the soil, 10 and 50 Percent, were considered. In Figure 49, where the soil’s moisture content is 10 percent, the output air temperature for the two-dimensional tube was higher by 0.84'F, than the . «...-pur— 87 .501.m:n. 0... .m c0_.m00..ux m> 00.3 on. c_ 05.300E3 :4 .wv 050.“. 2.8... 8:83-01 Dow 03.. mdm Dem mdm 0.6m m9 . 03 00 05.20:. o\oO— \ 05.205 $8 1 il 9 m— m... o. m. we: 5 3.8? .2 5%, £93 5 : om .0 83 892625-05 0m (21.) €9an 3H1 U! aJnleJedwai JIV 88 .83 9m 3.. m: 3 - -. m6; w> 93.32.25 I . L :< 50.30 . mg 059“. 3.30:. m8... mm . .. _ . 3 cm 9 3 2 o... No . . . . . . . . . . v.0 F.O . . . . . . m_ 82 oi u 3 82 A: - 9 m. (i) aJmeJadwei my 1nd1nQ 89 m.m .00 .. 03.. Q 0 0:0 0.. .0. 06; 0., 0.280080. :4 .3930 cm a .. . 05 .... $50.-.. 08; 3 mm 3 m: 2 3 S v.0 so 003. o; 003. Onw m— 3 m... @— 0566... 030m .0 m... 9 m. 003 c. >._00_0> 53. :05; 50:0. 9 .. On .0 003.. t (i) aJmeJadwei 19v 1nd1nQ 90 one-dimensional tube at 0.5 hour of the tube's operation. However, after 1 hour, the output air temperature for the two-dimensional tube was higher than the one-dimensional tube by l.3l°F. For the soil’s moisture content of 50 percent (see Figure 50), the output air temperature at 0.5 hour of the tube’s operation for two-dimensional tube, was higher than the one-dimensional tube by 1.59‘F. However, after 2 hours, the output air temperature for the two-dimensional tube was higher by 1.32°F, when compared to the one-dimensional tube. CHAPTER 5 CONCLUSIONS AND RECOMMENDATIONS 5.1 Conclusions Two mathematical models were considered, one a circular tube with radial coordinates and the other a square tube with Cartesian coordinates, known here as the one- and two-dimensional models. From the study that has been accomplished, several conclusions can be stated. Increasing the moisture content of the soil will increase the output air temperature of the tube. This is due to the increase in the heat capacity of the soil which depends on the moisture content of the soil. For the one-dimensional model with an air flow rate of 300 cubic feet per minute, when the moisture content in the soil is increased from 10 to 50 percent, the output air temperature of the tube is increased by 10.2 percent; with an air flow rate of 600 cubic feet per minute, the increment is 17.3 percent. As for the two-dimensional model, the increment is 0.2 percent when the air flow rate is 600 cubic feet per minute. It appears that with the increasing air flow rate, the increment gets larger. The formation of ice significantly influences the increase in temperature of the air in the tube. In the 91 one-d: per n perce hours tempe With outpu hours respe air 1 forma highs large stabi conte the ' tube outp‘ high hour the : tube aftEJ 92 one-dimensional tube with an air flow rate of 300 cubic feet per minute, when the moisture content of the soil is 10 percent, the output air temperature stabilized at about 12 hours of simulation time; and at 50 percent, the output air temperature continued to decrease gradually after 16 hours. With an air flow rate of 600 cubic feet per minute, the output air temperature stabilized at about 14 hours and 16 hours for 10 and 50 percent moisture content in the soil, respectively. This early stability phenomena in the output air temperature is due to the early occurance of the ice formation in the soil of 10 percent moisture content. For higher moisture content, additional time would be needed for larger ice formation and thus prolongs the occurance of the stabilized output air temperature. Therefore, the moisture content in the soil significantly affects the ice formation. It appears that the output air temperature for the two-dimensional tube is higher than the one-dimensional tube. For the 10 percent moisture content in the soil, the output air temperature for the two-dimensional tube is higher by 9.1 percent than the one-dimensional tube after 1 hour of operation. For the 50 percent moisture content in the soil, the output air temperature for the two-dimensional tube was higher than the one-dimensional tube by 9.2 percent after 2 hours of operation. 93 5.2 Recommendations The major problem encountered in this numerical study is that too much computer time is needed. This is due to the time step and space step being small. Therefore, some work should be done in modifying or refining the computer model. The program could be modified by having smaller space steps near the tube and larger ones away from the tube. This would reduce the amount of computer time without actually affecting the accuracy of the simulation. Another possibility to reduce the computer time is to use an approximate solution in calculating the ground temperature distribution. One area where the computer model could be refined is to include the tube's material and thickness, and the effects of the soil's thermal contact with the tube. Also, the calculation of the heat transfer coefficient of the soil to the atmospheric air at the ground surface should be included with real outside conditions of variable temperature and wind speed. Experimental research should be carried out for the validation of the numerical results to verify if the predicted results could actually be obtained. Experimental research should also be carried out to determine the soil's thermal properties since different references give different values of thermal properties for the same soil. The 8011’s thermal conductivity and heat capacity actually depend on its composition, moisture 94 content and temperature. Therefore, a suitable corelation should be obtained to include these effects. According to a parametric study (32) , the performance of the tube could be increased if the burial depth of the tube is increased. However, caution must be taken to avoid a permafrost region around the tube when in use in the northern climates. This numerical study is limited to the tube’s continuous operating condition. Study should be carried out in the "on/off" cyclic operating condition where air is allowed to enter the tube for 30 minutes and then shut off for the next 30 minutes for every hour. This is to allow the ground to recover between periods of the tube operation. Heat pumps are also used for summer cooling. Further study is called fOr, in simulating summer operation. For summer simulation, soil moisture migration effects should be included in the computer model because the moisture content in the soil plays an important role in the soil thermal properties. However, research in soil moisture migration with soil temperature has been limited. BIBLIOGRAPHY 10. BIBLIOGRAPHY Rice, C.K., Jackson, W.L., Fischer, S.K. and Ellison, R.D., "Design Optimization and the Limits of Steady-State Heating Efficiency for Conventional Single-Speed Air-Source Heat Pumps ", Oak Ridge, Tennessee: Oak Ridge National Laboratory, National Technical Information Service, 1981, Chapter 4. Gannon, Robert., "Ground-Water Heat Pumps", Popular Science, 1978. Same as 7(2) .‘ Bylinsky, Gene. , "Water to Burn", Fortune, October 20, 1980. Mei, V.C. and Fischer, S.I(., "A Theoretical and Experimental Analysis of Vertical, Concentric-Tube Ground-Coupled Heat Exchangers" , Oak Ridge, Tennessee: Oak Ridge National Laboratory, October, 1984. Parker, J .D. , Kavanaugh, S. and Ramanathan, R. , Performance Comparison of Air and Ground-Coupled Heat Pump Systems, Palo Alto, CaIifornia: Electric Power Research Institute,Final Report EM - 3408, 1984. Sumner, John A., Domestic Heat Pum s, Chalmington, Dorchester: Prism Press, SEaEeI Court. Dorset DT2 0H8, Great Britain, 1976. Metz, Philip 0., "Design, Operation and Performance of a Ground-Coupled Heat Pump System in a Cold Climate", Solar Technology Group, Brookhaven National Laboratory, Upton, New York, 1981. Mei, V.C. and Baxter, V.D., "Performance of a Ground-Coupled Heat Pump with Multiple Dissimilar U-Tube Coils in Series", Oak Ridge, Tennessee: Oak Ridge National laboratory. Rosenblad, 6., "Seasonal Heat Storing 1979-83 in Utby Ground Heat Pump Proj ect", International Conference 95 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. on Subsur_f;ace Heat Storage-Proceedings, Stockholm: Swedish Council for Bu ng Research, 1983, pp. 670-678. ' Mei, V.C., "Horizontal Ground-Coil Heat Exchanger Theoretical and Experimental Analysis", Oak Ridge, Tennessee: Oak Ridge Natinal Laboratory, December, 1986. Wasserman, David M,. "Evaluation and Simulation of a Crawl Space Coupled Heat Pump", M.S. Thesis, University of Tennessee, Knoxville, December, 1982. Robb, Ronald J ., "An Analytic Analysis of the Air-to-Air Ground-Coupled Heat Pump", M.S. Project Report, Michigan State University, East Lansing, March 9, 1983. Bose, J .E., "Earth Coil/Heat Pump Research at Oklahoma State University", Stillwater, Oklahoma: School of Technology, Oklahoma State University. Jumikis, Alfreds R., Thermal Geotechnics, New Brunswick, New Jersey: Rutgers University Press, 1977, Chapter 13. Thibault, Jules . , "Comparison of Nine Three-Dimensional Numerical Methods for the Solution of the ‘Heat Diffusion Equation", Numerical Heat Transfer, An International Journal of Computation and Methodology, 1985. Johnson, Lee W. and Riess, R. Dean. , Numerical Analysis, Addison-Wesley Publishing Company, Second E ition, 1982, Section 8.2.3. Johnson, Lee W. and Riess, R. Dean. , Nungrical Analysis, Addison-Wesley Publishing Company, Second Thibault, Jules. , "Comparison of Nine Three-Dimensional Numerical Methods for the Solution of the Heat Diffusion Equation", Numerical Heat Transfer, An International Journal of Computation and Methodology, 1985. Incropera, Frank P. and Dewitt, David P. , Introduction to Heat Transfer, John Wiley and Sons, New York, 1985, pp. 360. ' McAdams, W.H., Heat Transmission, Third Edition, Mcgraw-Hin, New York, 1954,4f1c apter 7. 96 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. Salmone, L.A., Kovacs, W.D. and Wichsler, H., ”Thermal Behavior of Fine-Grained Soils", Washington, D.C.: NBS Building Science Series - 149, 1982. Jennings, Burgess H., Environmental Engineering, Analysis and Practice, New York, New York: International TextboOk Company, 1970. U.S. Department of Commerce, National Oceanic and Atmospheric Administration , Environmental Data Service, in Cooperation with Michigan Weather Service. Same as (2). Koorevaar, P., Menelik, G. and Dirksen, C., Elements of Soil Physics, New York, New York: Elsevier Science Publishers, 1983. Johnson, Lee W. and, Riess, R. Dean., Numerical Analysis, AddisoneWesley Publishing Company, Second E ition, 1982, Chapter 8. Same as (23). Potter, Merle C. and Foss, John F., Fluid Mechanics, New York,New York: The Ronald Press Company,1975. Jumikis, Alfreds R., _Thermal Soil Mechanics, New Brunswick, New Jersey: Rutgers University Press, 1966, Chapter 2. Same as (11). Same as (11). 97 APPENDICES APPENDIX A COMPUTER.PROGRAM FOR ONE-DIMENSIONAL MODEL 000000000GOOOOOOOOOOOOOOOOOOOOOOOOO 98 Appendix A Computer Program for One-Dimensional Model PROGRAM ONEDIM *ttifitti************************************************************ This program is to determine the thermal performance of an under- ground heat tube in one dimension. It performs a finite difference analysis on a circular buried tube through which air is flowing. Central differences are used to calculate the heat gain of each element and forward differences are used in the time step. This program also calculates the thickness of ice formation around the tube. ******************************i*tiitt‘tit‘ki't**********************'k** Variable List TFAR TATM TEMP CPAIR H DNAIR AKSOLN AKSOLF CPSOLN CPSOLF CPSOIL DNSOIL CPFUSE CPWATR CPICE XINCRE RINCRE XNUM RNUM AMOIST VELCTY TIMLOP TIMEX Farfield Temperature of Soil Current Temperature of Atmosphere Array Containing Temperatures of Soil and Air Heat Capacity of Air Convective Heat Transfer Coefficient Density of Air Thermal Conductivity of Normal Soil Thermal Conductivity of Frozen Soil Heat Capacity of Soil in a Thawed State Heat Capacity of Soil in a Frozen State Heat Capacity of Soil based on its Current Temperature Density of Soil Heat of Fusion of Water Heat Capacity of Water Heat Capacity of Ice Increment Size in the x-direction Increment Size in the R-direction Number of Increments in the x-direction Number of Increments in the R-direction Moisture Content of Soil as a Percentage by Weight Velocity of Air Maximum Number of Timesteps to be Performed II 0000000000000000000000 00000 00000 99 AMASOL Dry Mass of Soil AMASAR Mass of Air DELI Array containing Change in Ice Formation RICE Array containing Ice Thickness RI Array containing Radius of Ice from Center of Tube AKMAJ Coefficient of Heat Transfer Equation AKMAJI II AKHAJZ ,, AKOO or AKI Coefficient of the J-1 term Heat Transfer Coefficient 2 Coefficient of the J+1 term Heat Transfer Coefficient Loop Counter Array Index II II WVH'UXQH; I! iii*********************i********fi**fiittiififiiiiii****i************** itiii*******i***iiiiiiififittifiiiii*iittiiiiittii********************* INTEGER I,J,K,T,P,XNUH,RNUH,TIMEX,TIHLOP,CHOICE,A,B,ANS,UNIT DIMENSION TEHP(12,10,2),R(10),RICE(11,2),DELI(11,2), +RI(11,2) ***************************t*i********i***i************************* Open Data File ii*********iifiiiitiittfiiiifiiiiifli*********************************** OPEN (UNIT'8,NAME"RESULT.DAT',TYPEB'UNKNOWN') tfi*************fi************iii*********i*************************** Initialize Test Variables iii*************ititiitiiiit*iittit**t****************************** 100 TFAR 39.4 TATM 0.0 CPAIR 0.24 DNAIR 0.08275 AKSOLN 0.453 ANSOLF 0.833 CPSOLN 0.3 CPSOLF 0.2 ONSOIL 95.0 CFFUSE 143.3 CPWATR 1.0 CFICE 0.5 PI 3.141593 K 1 XINCRE 5.0 XNUM 10 RNUM a R(2) 0.564 TIMLOP 12 TIMSTP 0.16667 TIMEX 12 C WRITE(*,*)' INPUT MOISTURE CONTENT OF THE SOIL' READ(*,*)AMOIST WRITE(8,*)' MOISTURE CONTENT OF SOIL - ',ANOIST WRITE(*,*)’ INPUT VELOCITY OF AIR IN THE TUBE’ READ(*,*)VELCTY WRITE(8,*)’ VELOCITY OF AIR IN THE TUBE - ',VELCTY WRITE(*,*)’ INPUT CONVECTIVE HEAT TRANSFER COEFFICIENT’. READ(*,*)H WRITE(8,*)’ CONVECTIVE HEAT TRANSFER COEFFICIENT = ’,H WRITE(*,*)' INPUT DELTA R’ READ(*,*)RINCRE WRITE(3,*)"OELTA R - ’,RINCRE C C *‘kfiii‘k‘ififiifiiifitfifiii*fliflfiiiiififiiiififiitiitflfii*itii’ifiitfiii'k'kitfiflukii**** 1131 C Set Initial Temperature c *iti*t***i******fl**fi*it**fifi***t**it***********it****************iiii C DO 100 I I 1,XNUM+2 DO 110 J I 1,RNUM+2 TEMP(I,J,K) I TFAR 110 CONTINUE 100 CONTINUE C c ***************************i**************************************** C Calculate Radius Of Nodes c ******************************ii*******i**************************** C WRITE(8,*)’RADIUS OE NOOES [R(3),R(4),.....]:’ DO 120 J I 3,RNUM+2 R(J) I R(J-l) + RINCRE WRITE(8,51)R(J) 51 FORMAT(1X,F8.5) 120 CONTINUE C c ****************i*****************t*ttiiii******t******************* C Calculate Time Step C ***********************************i***********t******************** C WRITE(8,54)TIMSTP 54 FORMAT(1X,’TIMSTP I ',F10.8) ******************i************************************************* Write Initial Ice Thickness, Radius of Ice and Temperature Of Air in the Tube ******************************************************************** O (I O C) O () WRITE(8,*)’ INITIAL ICE THICKNESS, RADIUS OP ICE AND’ WRITB(8,*)' TEMPERATURE OP AIR IN TEE TUBE :' DO 130 I - 2,XNUM+1 RICE(I,K) - 0.0 1132 RI(I,K) I R(2) DELI(I,K+1) I 0.0 WRITE(8,55)RICE(I,X),RI(I,K),TEMP(I,1,K) 55 FORMAT(1X,F10.8,3X,Flo.8,3X,F11.7) 130 CONTINUE C WRITB(8,2) 2 PORNAT(1x,70('*')) C C i*********t*****fiii*iititiitiiiittii*titttttttitit*i****tt*t******** c ******************************************************t************* c iii *** C *** Loop Calculating Heat Transfer Throughout The Array *** c it. *** C iii!it******t****ii*ifiti*fltittfliiiittittttii********i*************ia c tiiitititiiiiiifiiiititititiitfifi*itititttiti*fitiuittfltit**it**i*tit** C 00000 DO 200 T - 1,TINEx WRITE(8,*)’ T - ',T WRITE(8,3) EORNAT(1x,70('e')) DO 210 P - 1,TINLOP DO 230 I - 2,XNUN+1 DO 240 J - 2,RNUN+1 ***fiiiiiiifiitiiiififiiitfififiiiiiiit§fl**********i***t***i**************t Select Soil Heat Capacity ************i**fi******i***itttfifi******t*****************t*********** IF (TEMP(I,J,K) .GT. 32.0) THEN CPSOIL I CPSOLN + (ANOIST * CPWATR) ELSE 00000 00000 2103 IF (TEMP(I,J,K) .LT. 32.0) THEN CPSOIL I CPSOLF + (AMOIST * CPICE) ELSE CPSOIL I (CPSOLN+(AMOIST*CPWATR))*(RINCRE-RICE(I,K)) + /RINCRE + (CPSOLF+(AMOIST*CPICE))*RICE(I,K)/RINCRE END IF END IF *************i**t*i*fi**iifitifiiiiiiiii********i********************** Select Soil Thermal Conductivity *************************t**********tit*iit************************* TAVEl I (TEMP(I,J-1,K) + TENP(I,J,X)) / 2.0 TAVEZ I (TEMP(I,J+1,K) + TEMP(I,J,K)) / 2.0 IF (TAVEl .GE. 32.0) THEN AKSOLl I AKSOLN ELSE AKSOLI I AKSOLF END IF IE (TAVE2 .GE. 32.0) THEN AXSOLZ I AKSOLN ELSE AKSOLZ I AKSOLF END IF ********************************fi*********************************** Select Appropriate Heat Transfer Coefficient ******************************************************************** IE (J .E0. 2) THEN AMASOL - ONSOIL * PI * (R(2) + (RINCRE/4.0)) * RINCRE * + XINCRE AKMAJ - AMASOL * CPSOIL AKl - a * 4.0 * PI * R(2) . XINCRE AK2 - AKSOLZ . 4.0 * PI * (R(2) + (RINCRE/2.0)) * XINCRE / + RINCRE 00000 00000 104 CHOICE I 1 ELSE IF (J .20. 3) THEN AMASOL I DNSOIL * 2.0 * PI * R(3) * RINCRE * XINCRE AKMAJ I AHASOL * CPSOIL AKl - AKSOL1 . 2.0 a PI * (R(3) - (RINCRE / 2.0)) * XINCRE/ + RINCRE AKZII.AKSOL2 * 2.0 * PI * (R(3) + (RINCRE / 2.0)) * XINCRE/ + RINCRE CHOICE I 2 ELSE ANASOL - ONSOIL * 2.0 . PI . R(J) . RINCRE * XINCRE AKHAJ - AMASOL . CPSOIL AKl - AKSOLI * 2.0 * PI * (R(J) - (RINCRE / 2.0)) . xINCRE/ + RINCRE AKZ - AKSOL2 * 2.0 9 PI . (R(J) + (RINCRE / 2.0)) * XINCRE/ + RINCRE CHOICE - 3 END IF ********fi********i******tiiiiifitfiitiiifiiiiit*******i*i************** Calculate The Temperature of The Element *it*********t***fiititttfitttfit*ttt*****t***************************** -TEMP(I,J,K+1) I TIMSTP * ( AKl * (TEMP(I,J-1,K) - TEMP(I,J,K)) + ‘ + AKZ * (TEMP(I,J+1,K) - TEMP(I,J,K))) + / mm + TENP(I,J.K) **********i********i**tit******************************************* Calculating Ice Formation and Making Temperature Adjustments **********************************ti***iittit**********************t IF (CHOICE .E0. 1) THEN IP (TEMP(I,J,K+1) .LT. 32.0) THEN IE (RICE(I,K) .LE. (RINCRE/2.0)) THEN 105 AKMAJZ I DNSOIL * CPFUSE * AMOIST * 2.0 * + PI * RI(I,K) * XINCRE AKOO I AHASOL * CPSOIL / TIMSTP DELI(I,H+1) - TIHSTP * + (ARGO * (TEMP(I,J,K+1) - TEMP(I,J,K))) + . / AKHAJZ IF (DELI(I,K+1) .LT. 0.0) THEN DELI(I,K+1) I -DELI(I,K+1) RICE(I,K+1) I RICE(I,K) + DELI(I,K+1) RI(I,K+1) I 3(2) + RICE(I,K+1) ELSE RICE(I,K+1) I RICE(I,K) + DELI(I,K+1) RI(I,K+1) I R(2) + RICE(I,K+1) END IF IE (RICE(I,K+1) .LT. (RINCRE/2.0)) THEN TEMP(I,J,H+1) . 32.0 END IE END IE ELSE RICE(I,K+1) - RICE(I,K) RI(I,K+1) - R(2) + RICE(I,K+1) END IE ELSE IE (CHOICE .Eq. 2) THEN IE (RICE(I,K) .GT. (RINCRE/2.0)) THEN AKHAJZ - DNSOIL . CPEUSE . ANOIST . 2.0 * + PI * RI(I,K) * XINCRE AKOO - AHASOL * CPSOIL / TIMSTP DELI(I,K+1) - TINSTP * (Axoo * (TEMP(I,J,K+1) - TEMP(I,J,K))) / ARMAJz IF(DELI(I,K+1) .LT.0.0)THEN 1136 DELI(I,K+1)I-DELI(I,K+1) RICB(I,K+1) - RIC£(I,K) + DELI(I,K+1) RI(I,H+1) - R(2) + RICE(I,R+1) END IE . . - r v -. ‘- I. - .- J. v v - e- a .... ». r‘. ... . — "I. A .a - s .- -- IE (RICE(I,H+1) .LT. (3.0*RINCRE/2.0)) THEN IE (RICE(I,K+1) .GT. (RINCRE/2.0)) THEN IE (TEHP(I,J,K+1) .LT. 32.0) THEN TENP(I,J,H+1) - 32.0 4 END IE ‘ .' ' END IF END IF END IF END IF 240 CONTINUE C 230 CONTINUE *************i************fi****fi**t********************************* Reset Tube Entrance Temperature to Current Ambient Temperature iit*iit******ttiiii**i**t*i****i******************fi***************tt 00000 TEMP(1,1,K+1) I TATM i*******************i*********************t************************* Match End Wall Temperature to Satisfy the Insulation Condition ttit*t*************************i**fi********************************* 00000 TEMP(1,2,K+1) I TEMP(2,2,K+1) 0 C ****************i**fi****fifi*i*fi*****fi*******fi**fi********************* 107 C Calculate The Temperature of Air in The Tube C ****i*******ifi**itttfiiifliiifitfitttfitfiiiitfifi**i*t********************* C AKMAJI I (2.0 * H * XINCRE) / (DNAIR * CPAIR * R(2) * VELCTY) DO 220 I I 2,XNUH+1 TEMP(I,I,K+1)IAKMAJ1 * (TEMP(I-l,2,K+l) - TEHP(I-1,1,x+1)) + + TEHP(I-1,1,K+l) 220 CONTINUE C C ******************************************************************** C Move K+1 Temperatures to K for Next Loop C ************i****iifiifiii*****************i************************** C 00 310 I I 2,XNUH+1 DO 320 J I 1,RNUH+1 TEMP(I,J,K) I TmP(I,J,K+l) 320 CONTINUE 310 CONTINUE C C **i**************fi********i***i*********i*************************** C Move K+1 Ice Formation to K for Next Loop C *tfii**********ttfi**iiflitiii*tfltflfifltiiifitifliittfi*itii*ititii*****i*** C DO 330 I I 2,XNUM+1 RICE(I,X) I RICE(I,X+1) RI(I,K) . RI(I,H+1) 330 CONTINUE C C 210 CONTINUE C C *i*********i*fi*t****itiiitt*fitiiitfiiiitttiiiiifiittiit*************** C Write Out The Temperature of Soil Elements C *******i****i******************************************************* C WRITE(8,*)’ TEMPERATURE OF SOIL ELEMENTS FOR EACH X-INCREMENT: 2108 +1 DO 500 J I 2,RNUM+2 WRITE(8,13)TEMP(2,J,K),TEMP(3,J,K),TEMP(4,J,K),TEMP(S +.J.K).TEMP(6.J.K) 13 EORMAT(1x,5E11.7) 500 CONTINUE WRITE(8,*) C DO 501 J I 2,RNUM+2 WRITE(8,13)TEMP(7,J,K),TEMP(8,J,K),TEMP(9,J,K), +TEMP(10,J,K),TEMP(11,J,K) 501 CONTINUE WRITE(8,*) C C **************iiiiiitiiitfiiiti*§§****fi****************************** C Write Out The TemperatUre of Air in The Tube C it*****************tttittflttitiitfittitti*********t****************** C WRITE(8,*)' TEMPERATURE OF AIR IN THE TUBE: ' Do 225 I I 2.6 WRITE(8,14)TEMP(I,1,K+1),TEMP(I+5,1,K+1) 14 FORMAT(1X,2F11.7) 225 CONTINUE C C ***********************************t******************************** C Calculate Time for Every TIMLOP C ******************i**t****i***************************************** C HOURS I TIMSTP * TIMLOP * T WRITE(*,16)TEMP(XNUM+1,1,K+l) WRITE(8,16)TEMP(XNUM+1,1,K+1) 16 EORMAT (1x,' OUTPUT TEMPERATURE OE TUBE Is ’,Fll.7) WRITE(*,17)HOURS WRITE(8,17)HOURS 17 FORMAT(1X,' TIME Is ’,F12.8,’ HOURS’) 109 C **********fi********fi**i****§*ifl**ittitfliflttii**********************t C Write Out Ice Formation For Every TIMLOP C *********it*****i*tti***tttttittittiitiittit*fi********************** WRITE(8,*) WRITE(8,*)’ THICKNESS OF ICE FORMED:' DO 600 I I 2.6 WRITE(8,*)RICE(I,K),RICE(I+5,K) 600 CONTINUE C C 200 CONTINUE C C *****************tt************************************************* C If Maximum Time Step is Reached Write Out Message C t**********ititifiiiittitfitiitttttittfiitiiifititiititiitt************* C WRITE(8,18) 18 FORMAT(1X,70('*’)) NRITE(8,19) 19 FORMAT(1X,70(’*')) WRITE(*,*)’ MAXIMUM TIME EXCEEDED' WRITE(8,*)' MAXIMUM TIME EXCEEDED’ C C ii***********i************§*********************tit*i*************** C If Stable Temperature is Achieved, Write Out Message C ii*iititiifiiiiifiiifliifiii******************************************** C 420 WRITE(*,*)’ PROGRAM TERMINATED’ WRITE(8,*)' PROGRAM TERMINATED' CLOSE(8) STOP END APPENDIX B COMPUTER PROGRAM FOR TWO-DIMENSIONAL MODEL 00000000000000000000000000000000000 111) Appendix B Computer Program for Two-Dimensional Model PROGRAM TWODIM iiiiii*iiiflflfifiififlifiifififififiifififlfiit...iii...*****i****fi**************** This program is to determine the thermal performance of an under ground heat tube in two dimension. It performs a finite difference analysis on a square buried tube through which air is flowing. Central differences are used to calculate the heat gain of each element and forward differences are used in the time step. The program also calculates the thickness of ice formation around the tube. ******fi***********i*********fi************i************************** variable list TSUR TATM TSOILD TEMPT TEMP TDIFF TDIFES CPAIR DNAIR AKSOLN AKSOLF CPSOLN CPSOLE CPSOLC DNSOIL CPFUSE CPWATR CPICE XINCRE YINCRE ZINCRE XNUM YNUM soil surface temperature current temperature of the atmosphere deep soil temperature temporary temperature of preprocessed increments array containing temperatures of soil and air temperature difference temperature difference required for stability heat capacity of air density of air thermal conductivity of normal soil thermal conductivity of frozen soil heat capacity of soil in a thawed state heat capacity of soil in a frozen state heat capacity Of soil in a freezing state density Of soil heat of fusion Of water heat capacity of water heat capacity of ice the increment size in the x-direction the increment size in the y-direction the increment size in the z-direction number of increments in the x-direction number of increments in the y-direction 000000000000000000000000000000000000 m >'ra m t* x C4 H ZNUM TNUM VELCTY AMOIST HSUR TIMEX TIMLOP AMASOL AMASAR TUBEHT DELIl DELIZ HICE1 HICEZ 111. number of increments in the z-direction the depth at which the tube is located vertically the velocity of air moisture content Of soil as a percentage by weight soil-air heat transfer coefficient at the surface soil-air heat transfer coefficient at the tube maximum number of time steps to be performed II mass of dry soil mass of air location of tube vertically array containing change in ice formation above the tube array containing change in ice formation below the tube array containing ice formation above the tube array containing ice formation below the tube coefficient of heat transfer equation 0' II coefficient of the J-1 term coefficient of the J+l term coefficient of the K-1 term coefficient of the x+1 term coefficient of the air term II energy released/absorbed due to formation/defrosting of ice lOOp counter array index ' II II II II II II II fl 00000 00000 00000 2112 ********************i*fl*tfifiifiifiittfiifififii*iitiifi***************i***** ************************fi**********.****fl**************************i INTEGER I,J,K,L,T,P,A,B,C,XNUM,YNUM,ZNUM,TNUM,TUBEHT,TIMEX, +TIMLOP,CHOICE,ANS,UNIT ' DIMENSION TEMP(52,17,26,2),HICE1(51,2),HICE2(51,2),DEL11(51,2) +,DELI2(51,2) *ii*********i****i*i**i**************i**********t******************* open data file *i******************t*******t*****i*******i*************t*********** OPEN (UNITI8,FILEI’RESULT.DAT',STATUSI’NEW') ***********i************§*i*§**i*****i****************************** initialize test variables ***i**********iii*******Qt***i**ttfii*tfifiiitiit*********************i TSUR I 30.0 TATM I 0.0 TSOILD I 39.7 CPAIR I 0.24 DNAIR I 0.08275 AKSOLN I 0.458 AKSOLF I 0.833 CPSOLN I 0.3 CPSOLF I 0.2 DNSOIL I 95.0 CPFUSE I 143.3 CPWATR I 1.0 CPICE I 0.5 XINCRE I 1.0 YINCRE I 0.5 ZINCRE I 0.5 000000 120 110 .113 XNUM I 50 YNUM I 14 ZNUM I 24 TNUM I 10 L I 1 HSUR I 0.959 WRITE(*,*)’ INPUT MOISTURE CONTENT OE SOIL' READ(*,*)AMOIST WRITE(8,*)’ MOISTURE CONTENT OE SOIL - ',AMOIST WRITE(*,*)’ INPUT TEMPERATURE DIEEERENCE EOR STABILITY' READ(*,*)TDIFFS WRITE(8,*)’ TEMPERATURE DIEEERENCE EOR STABILITY - ’,TDIFFS WRITE(*,*)' INPUT VELOCITY OE AIR IN THE TUBE' READ(*,*)VELCTY WRITE(8,*)’ VELOCITY OE AIR IN THE TUBE - ’,VELCTY WRITE(*,*)' INPUT CONVECTIVE HEAT TRANSPER COEFFICIENT’ READ(*,*)H WRITE(8,*)’ CONVECTIVE HEAT TRANSEER COEEEICIENT - ’,H WRITE(*,*)’ INPUT TIMLOP’ READ(*,*)TIMLOP WRITE(8,*)’ TIMLOP - ',TIMLOP ***************************t**************************************** set temperature of soil to its farfield temperature - the air temperature in the tube is set to soil temperature at that depth ******************fi***i***********t********************************* DO 100 K I 1,2NUM-2 TEMPT I TSOILD - ((TSOILD - 32.0) * (0.001959 * FLOAT(K) +**2.0 - 0.001558 * FLOAT(K) - 0.000401)) DO 110 I I 1,XNUM+2 DO 120 J I 1,YNUM+3 TEMP(I,J,K,L) I TEMPT CONTINUE CONTINUE 100 150 140 130 170 160 70 80 185 114 CONTINUE DO 130 K I ZNUM-I,ZNUM+1 TEMPT - (TSUR - 32.0)/2.0 . ELOAT(H) - (23.0 / 2.0 . TSUR) ++400.0 DO 140 I I l,XNUM+2 DO 150 J I l,YNUM+3 TEMP(I,J,X,L) I TEMPT CONTINUE CONTINUE CONTINUE K I ZNUM+2 DO 160 I I 1,XNUM+2 DO 170 J I 1,YNUM+3 TEMP(I,J,K,L) I TATM CONTINUE CONTINUE WRITE(8,*) WRITE(8,*)’ INITIAL TEMPERATURE OE THE NODEs POR ALL x :’ WRITE(8,15) DO 180 I - 1,52 WRITE(8,*) NRITE(8,70) FORMAT(1X,'K',4X,’J - l',5X,’J - 2’,5X,’J - 3’,5X,’J = 4', +5X,’J - 5’,5X,'J - 6’) DO 135 K - zNUM+2,1,-1 ‘ WRITE(8,80)K,TEMP(I,1,K,L),TEMP(I,2,K,L),TEMP(I,3,K,L), +TEMP(I,4,X,L),TEMP(I,5,K,L),TEMP(I,6,K,L) EORMAT(1x,Iz,SE11.7) CONTINUE WRITE(8,*) WRITE(8,71) 2115 71 FORMAT(1X,’K',4X,'J I 7’,5X,’J I 8’,SX,’J I 9’,4X,’J I 10’, +4X,’J I 11’,4X,’J I 12’) DO 190 K I ZNUM+2, 1,-1 WRITE(8,80)K,TEMP(I,7,K,L),TEMP(I,8,K,L),TEMP(I,9,K,L), +TEMP(I,10,K,L),TEMP(I,11,K,L),TEMP(I,12,K,L) 190 CONTINUE C WRITE(8,*) WRITE(8,72) 72 FORMAT(1X,’X’,3X,’J I 13’,4X,’J I 14’,4X,’J I 15’,4X, +’J - 16’,4X,’J - 17') DO 195 K - zNUM+2,1,-1 WRITE(8,81)X,TEMP(I,13,K,L),TEMP(I,l4,K,L), +TEMP(I,15,H,L),TEMP(I,16,H,L),TEMP(I,17,H,L) 81 EORMAT(1x,I2,5E11.7) 195 CONTINUE C 180 CONTINUE ***********t*fi********fifii*fi*i******iii****************************** calculate heat capacity of freezing soil ***i********flfifi**t*iiflififlfltiiiiiifiiiii*********ifi******************* 00000 CPSOLC I CPSOLN + (AMOIST * CPWATR) + (AMOIST * CPFUSE / 2.0) WRITE(8,3)CPSOLC FORMAT(1X,’ HEAT CAPACITY OF FREEZING SOIL I ’,F7.3) ****************i****i*****i***********************************i**** calculate time step **************i***************************************************** 00000“ TIMSTP - XINCRE / VELCTY WRITE(8,4)TIMSTP 4 FORMAT(1X,' TIME STEP . ’,F10.8) C *ifli***************iiifi*****t****fifliiti**t***fi*fi******************** 3116 C calculate tube location underground C ******************t***t**t*****ttit*itittttttt********************** C TUBEHT I ZNUM - TNUM + 1 C C ******ii*iiiiifii‘h*****fl**fi**fi*i***¢**iiii“.**************i*********** C calculate mass of air in the tube C **ti************i***ifiii******************************************** C AMASAR I 4.0 * DNAIR * XINCRE * YINCRE * ZINCRE C C *****i************************************************************** C calculate mass of soil in an element C ****************i******************ii******t************************ C AMASOL I DNSOIL * XINCRE * YINCRE * ZINCRE C C ***********************i*******i************************************ C write on-data file initial temperature of air along the tube C ********i***********tfiitiitiii*********************i***********i**** C WRITE(8,*) WRITE(8,*)’ INITIAL TEMPERATURE OF AIR ALONG THE TUBE :’ WRITE(8,15) WRITE(8,5) 5 FORMAT(1X,’I’,1X,'TEMPERATURE’) DO 200 I I 2,XNUM+1 WRITE(8,6)I,TEMP(I,2,TUBEHT,L) 6 FORMAT(1X,IZ,1X,F11.7) 200 CONTINUE 00000 ******************************************************************** set ice formation around the tube to be zero initially **************tiiifiii**********iiit********************************* .117 WRITE(8,*) WRITE(8,*)’ INITIAL ICE THICKNESS :' WRITE(8,15) . WRITE(8,7) 7 FORMAT(1X,’I’,4X,'HICEl',SX,'HICEZ’) DO 210 I - 2,XNUM+1 HICE1(I,L) - 0.0 HICE2(I,L) - 0.0 HICE1(I,L+1) - 0.0 HICE2(I,L+1) - 0.0 DELIl(I,L+l) - 0.0 DELI2(I,L+1) - 0.0 WRITE(8,8)I,HICE1(I,L),HICE2(I,L) 8 EORMAT(1x,I2,2E10.B) 210 CONTINUE C C WRITE(8,10) WRITE(8,lO) 10 FORMAT(1X,70('*’)) C C titttfittitttt*titttttitttttttttttiitiitiittttittitaitiii C *********iiiiiitttttfiitttitittt******************************* c ****** ******* C ****** LOOP CALCULATING HEAT TRANSEER THROUGHOUT THE ARRAY ******* c ****** ******* C i********************itiiitttittttttt*tt'k'i******************** C *ititfitiiitttittittitttttttt*ittiittttttittitit********* C 20 WRITE(*,*)’ INPUT TIMEX’ READ(*,*)TIMEX WRITE(8,*)' TIMEx - ’,TIMEX DO 300 T - 1,TIMEx WRITE(8,*)’ T - ’,T WRITE(8,15) 15 EORMAT(1x,70('-’)) f 1118 DO 310 P I 1,TIMLOP DO 320 I I 2,XNUM+1 DO 330 J I 2,!NUM+2 DO 340 K I 2,ZNUM+1 select soil heat capacity C C c ****************t*titfitifittittttiitit**********************i******** C c t***********titttttttttittttitititiiititiifittii*******************i* IF (TEMP(I,J,K,L) .GT. 32.0) THEN CPSOIL I CPSOLN + (AMOIST * CPWATR) ELSE IF (TEMP(I,J,K,L) .DT. 32.0) THEN CPSOIL I CPSOLF + (AMOIST * CPICE) ELSE CPSOIL I CPSOLC END IF END IF C C ********i***fiiflfiiifiitifliiiiiiiifiiiii******************f************* C select soil thermal conductivity C iti**********i*iiiitiitifittitititiiit******************************* C TAVEl I (TEMP(I,J-1,K,L) + TEMP(I,J,K,L)) / 2.0 TAVE2 I (TEMP(I,J+1,K,L) + TEMP(I,J,K,L)) / 2.0 TAVE3 I (TEMP(I,J,K-1,L) + TEMP(I,J,R,L)) / 2.0 TAVE4 I (TEMP(I,J,K+1,L) + TEMP(I,J,K,L)) / 2.0 C IF (TAVEl .GE. 32.0) THEN AKSOLl I AKSOLN ELSE AKSOLI I AKSOLF END IF C IF (TAVE2 .GE. 32.0) THEN AKSOLZ I AKSOLN 2119 ELSE AKSOLZ I AKSOLF END IF C ' IE (TAVE3 .GE. 32.0) THEN AKSOL3 I AKSOLN ELSE AKSOLJ I AKSOLP END IF C IF (TAVE4 .GE. 32.0) THEN AKSOLA I AKSOLN ELSE AKSOL4 I AKSOLF END IF C C *iit*i**************i**t************************'k C ************i*t******iiiit**********************ttitittti'ki'kt******* C ** select the appropriate heat transfer coefficients ** C ******i*******tii**iti*ttiitttfiiiiiiitt*****i*********************** C **************t********************************** C .C if node borders the center of bottom face of tube C IF (J.EQ.2 .AND. K.EQ.TUBEHT-l) THEN AKMAJ I TIMSTP / (AMASOL * CPSOIL) AK1 I AKSOLI * XINCRE * ZINCRE / YINCRE AKZ I AKSOLZ * XINCRE * ZINCRE / YINCRE AKB I 2.0 * AKSOLJ * XINCRE * YINCRE / ZINCRE AK4 I 2.0 * H * XINCRE * YINCRE AKS I 0.0 AK6 I 0.0 DELTAE I 0.0 CHOICE I l C C if node is directly below the node which borders the center of 12C) C bottom face of tube C ELSE IF (J.EQ.2 .AND. K.EQ.TUBEHT-2) THEN AKMAJ I TIMSTP / (AMASOL * CPSOIL) AKI I AKSOLI * XINCRE * ZINCRE / YINCRE AK2 I AKSOL2 * XINCRE * ZINCRE / YINCRE A83 I AKSOL3'* XINCRE I YINCRE / ZINCRE AK4 I AKSOL4 * XINCRE * YINCRE / ZINCRE AKS I 0.0 AK6 I 0.0 DELTAE I 0.0 CHOICE I 2 C C if node borders the center of top face of tube C ELSE_IF (J.EQ.2 .AND. K.EQ.TUBEHT+1) THEN AKMAJ I TIMSTP / (AMASOL * CPSOIL) AMI I AKSOLl * XINCRE * ZINCRE / YINCRE AKZ I AKSOL2 * RINCRE * ZINCRE / YINCRE AK3 I 2.0 * H * XINCRE * YINCRE AK4 I 2.0 * AKSOL4 * XINCRE * YINCRE / ZINCRE AKS I 0.0 AKG I 0.0 DELTAE I 0.0 CHOICE I 3 C C if node is directly above the node which borders the center of C top face of tube C ELSE IF (J.EQ.2 .AND. K.EQ.TUBEHT+2) THEN AKMAJ I TIMSTP / (AMASOL * CPSOIL) AKl I AKSOLI * XINCRE AH2 I AKSOLZ * XINCRE AK3 I AKSOLJ * XINCRE AK4 I AKSOL4 * XINCRE AKS I 0.0 * ZINCRE / YINCRE * ZINCRE / YINCRE * YINCRE / ZINCRE * YINCRE / ZINCRE AK6 I 0.0 DELTAE I 0.0 CHOICE I 4 C 1121 C if node borders the bottom corner of tube C ELSE IF (J.EQ.3 IE (HICE2(I,L) .LT. (0.5*ZINCRE)) THEN C there could be some freezing, defrosting or no freezing at all AKMAJ I 2.0 * TIMSTP / (3.0 * AMASOL * CPSOIL) AKl I AKZ I AK3 I AK4 I AKS I AKG I 0.0 0.0 0.0 0.0 0.0 0.0 .AND. X.EQ.TUBEHT-1) THEN DELTAE I DNSOIL I CPFUSE * AMOIST * XINCRE * +DELIZ(I,L+1) * (YINCRE + ZINCRE + 4.0 * HICE2(I,L) + +2.0 * DELI2(I,L+1)) / TIMSTP CHOICE I 5 ELSE IF (HICEZ(I,L) .EQ. (0.5*ZINCRE)) THEN IF (HICE2(I,L+l) .LT. (O.5*ZINCRE)) THEN C some defrosting occurs I 2.0 * TIMSTP / (3.0 * AMASOL * CPSOIL) 0.0 0.0 0.0 0.0 0.0 0.0 DELTAE I DNSOIL * CPFUSE * AMOIST * XINCRE * +DELI2(I,L+1) * (YINCRE + ZINCRE + 4.0 * HICE2(I,L) + +2.0 * DELIZ(I,L+1)) / TIMSTP CHOICE I 6 2122 ELSE C element is completely frozen C element is completely frozen 0 ()(3 0 AKMAJ - 2.0 . TIMSTP / (3.0 * AMASOL . CPSOIL) ARI - ARSOLI . XINCRE * zINCRE / YINCRE sz - 2.0 . ARSOL2 . XINCRE . ZINCRE / YINCRE AH3 - 2.0 . AHSOL3 « xINCRE * YINCRE / zINCRE AK4 - ARSOL4 . xINCRE . YINCRE / zINCRE - AKS H a XINCRE * zINCRE AK6 - H . XINCRE . YINCRE DELTAE - 0.0 CHOICE - 7 ‘~ END IE ELSE AKMAJ I 2.0 * TIMSTP / (3.0 * AMASOL * CPSOIL) ARI I AKSOLI * XINCRE * ZINCRE / YINCRE AKZ I 2.0 * AKSOLZ * XINCRE * ZINCRE / YINCRE AK3 I 2.0 * AKSOL3 * XINCRE * YINCRE / ZINCRE AK4 I AKSOL4 * XINCRE * YINCRE / ZINCRE AKS I H * XINCRE * ZINCRE AKS I H * XINCRE * YINCRE DELTAE I 0.0 CHOICE I 8 END IF if node is directly below the node which borders the bottom corner of tube ELSE IF (J.EQ.3 .AND. K.EQ.TUBEHT-2) THEN IE (HICE2(I,L) .LE. (O.5*ZINCRE)) THEN IE (HICE2(I,L+1) .GT. (O.S*ZINCRE)) THEN there is some freezing occuring AKMAJ I TIMSTP / (AMASOL * CPSOIL) AKl I 0.0 AKZ I 0.0 123 AH3 - 0.0 AK4 - 0.0 Axs - 0.0 AK6 - 0.0 DELTAE - DNSOIL . CPEUSE * AMOIST . XINCRE * YINCRE * +OELIz(I,L+1) / TIMSTP CHOICE - 9 c ELSE C no freezing occurs ' AKMAJ - TIMSTP / (AMASOL * CPSOIL) AKl - AHSOLI . xINCRE . ZINCRE / YINCRE sz - AKSOLZ . XINCRE * ZINCRE / YINCRE Ax: - AHSOLS . XINCRE * YINCRE / ZINCRE AK4 - AKSOL4 . XINCRE . YINCRE / ZINCRE AKS - 0.0 AKG - 0.0 DELTAE - 0.0 CHOICE . 10 END IE ELSE C there could be some freezing or defrosting AKMAJ I TIMSTP / (AMASOL * CPSOIL) AKl I 0.0 AK2 I 0.0 AK3 I 0.0 AK4 I 0.0 AKS I 0.0 AK6 I 0.0 DELTAE I DNSOIL*CPFUSE*AMOIST*XINCRE*YINCRE* +DELIZ(I,L+l)/TIMSTP CHOICE - 11 END IF C C if node borders the center of side face of tube 2124 ELSE IF (J.EQ.3 .AND. K.EQ.TUBEHT) THEN IF (HICEl(I,L) .LT. (O.5*ZINCRE)) THEN C there could be some freezing, defrosting or no freezing at all AKMAJ I TIMSTP / (AMASOL * CPSOIL) AKl I 2.0 I H * XINCRE * ZINCRE AKZ I 2.0 * AKSOLZ I XINCRE * ZINCRE / YINCRE AK: I AKSOLJ * XINCRE * YINCRE / ZINCRE AK4 I AKSOL4 * XINCRE * YINCRE / ZINCRE AKS I 0.0 AKG I 0.0 DELTAE I DNSOIL * CPFUSE * AMOIST * XINCRE * +SQRT(ZINCRE ** 2.0 + ABS(HICEI(I,L) - HICE2(I,L)) ** 2.0) * +(DELII(I,L+1) + DELIZ(I,L+1)) / TIMSTP CHOICE I 12 ELSE IF (HICE1(I,L) .EQ. (0.S*ZINCRE)) THEN IF (HICE1(I,L+1) .LT. (0.5*ZINCRE)) THEN C some defrosting occurs AKMAJ I TIMSTP / (AMASOL * CPSOIL) ARI I 0.0 AK2 I 0.0 AK3 I 0.0 AK4 I 0.0 AKS I 0.0 AK6 I 0.0 DELTAE I DNSOIL * CPFUSE * AMOIST * XINCRE * +SQRT(ZINCRE M 2.0 + ABS(HICE1(I,L) - HICE2(I,L)) u 2.0) * +(DELI1(I,L+1) + DELIZ(I,L+1)) / TIMSTP CHOICE I 13 ELSE C element is completely frozen AKMAJ I TIMSTP / (AMASOL * CPSOIL) AKl I 2.0 * H * XINCRE * ZINCRE AKZ I 2.0 * AKSOLZ * XINCRE * ZINCRE / YINCRE 1255 AK3 I AHSOL3 * XINCRE * YINCRE / ZINCRE AK4 I AKSOL4 I XINCRE * YINCRE / ZINCRE AKS I 0.0 AKG I 0.0 DEDTAE I 0.0 CHOICE I 14 END IF {- ELSE C element is completely frozen AKMAJ I TIMSTP / (AMASOL * CPSOIL) AKl I 2.0 * H * XINCRE * ZINCRE AKZ I 2.0 * AKSOLZ * XINCRE * ZINCRE / YINCRE AK: I AKSOLJ I XINCRE * YINCRE / ZINCRE AK4 I AKSOL4 * XINCRE I YINCRE / ZINCRE AKS I 0.0 AKS I 0.0 DELTAE I 0.0 CHOICE I 15 END IF C C if node borders the top corner of tube C ELSE IF (J.EQ.3 .AND. K.EQ.TUBEHT+1) THEN IF (HICE1(I,L) .LT. (0.5*ZINCRE)) THEN C there could be some freezing, defrosting or no freezing at all AKMAJ I 2.0 * TIMSTP / (3.0 * AMASOL * CPSOIL) AKI I AKSOLl * XINCRE * ZINCRE / YINCRE AKZ I 2.0 * AKSOLZ * XINCRE * ZINCRE / YINCRE AKB I AKSOLB * XINCRE * YINCRE / ZINCRE AK4 I 2.0 * AKSOL4 * XINCRE * YINCRE / ZINCRE AKS I H * XINCRE * ZINCRE AKG I H * XINCRE * YINCRE DELTAE I DNSOIL * CPFUSE * AMOIST * XINCRE * +DELIl(I,L+1) * (YINCRE + ZINCRE + 4.0 * HICE1(I,L) + +2.0 * DELIl(I,L+l)) / TIMSTP CHOICE I 16 1265 ELSE IF (HICE1(I,L) IF (HICE1(I,L+1) C some defrosting occurs I 2.0 * TIMSTP / (3.0 * AMASOL * CPSOIL) AKMAJ AKl I 0.0 AX2 I 0.0 AK3 I 0.0 AK4 I 0.0 AKS I 0.0 AKG I 0.0 .EQ. (0.5*ZINCRE)) THEN .LT. (0.5*ZINCRE)) THEN DELTAE - DNSOIL . CPEUSE . AMOIST * XINCRE * +DELI1(I,L+1) . (YINCRE + zINCRE + 4.0 . HICE1(I,L) + +2.0 * DELIl(I,L+l)) / TIMSTP CHOICE - 17 ELSE C element is completely frozen AKMAJ AKG END IF ELSE I 2.0 * TIMSTP / (3.0 * AMASOL * CPSOIL) ARI I AKSOLl * XINCRE * ZINCRE / YINCRE 2.0 * AKSOLZ * XINCRE * ZINCRE / YINCRE AKSOLJ * XINCRE * YINCRE / ZINCRE I 2.0 * AXSOL4 * XINCRE * YINCRE / ZINCRE I H * XINCRE * ZINCRE H * XINCRE * YINCRE DELTAE I 0.0 CHOICE I 18 I C element is completely frozen AKMAJ I 2.0 * TIMSTP / (3.0 * AMASOL * CPSOIL) ARI I AKSOLl * XINCRE * ZINCRE / YINCRE AKZ I 2.0 * AKSOL2 * XINCRE * ZINCRE / YINCRE AKJ I AKSOL3 * XINCRE * YINCRE / ZINCRE 127' AK4 I 2.0 * AKSOL4 * XINCRE * YINCRE / ZINCRE AKS I H * XINCRE * ZINCRE ANS I H * XINCRE * YINCRE DELTAE I 0.0 CHOICE I 19 END IF if node is directly above the node which borders the top corner of tube 0000 ELSE IF (J.EQ.3 .AND. K.EQ.TUBEHT+2) THEN IF (HICE1(I,L) .LE. (0.5*ZINCRE)) THEN IF (HICE1(I,L+1) .GT. (0.5*ZINCRE)) THEN C there is some freezing occuring AKMAJ I TIMSTP / (AMASOL * CPSOIL) AKl I 0.0 AK2 I 0.0 AK3 I 0.0 AK4 I 0.0 AKS I 0.0 AKS I 0.0 DELTAE I DNSOIL * CPFUSE * AMOIST * XINCRE * YINCRE * +DELIl(I,L+l) / TIMSTP CHOICE I 20 ELSE C no freezing occurs AKMAJ I TIMSTP / (AMASOL * CPSOIL) AKl I AKSOLl * XINCRE * ZINCRE / YINCRE AKZ I AKSOLZ * XINCRE * ZINCRE / YINCRE AK3 I AKSOL3 * XINCRE * YINCRE / ZINCRE AK4 I AKSOL4 * XINCRE * YINCRE / ZINCRE AKS I 0.0 AK6 I 0.0 DELTAE I 0.0 CHOICE I 21 128 END IF C ELSE C there could be some freezing or defrosting AKMAJ I TIMSTP / (AMASOL * CPSOIL) AH1 - 0.0 j AKZ I 0.0 AK3 I 0.0 AK4 I 0.0 AKS - 0.0 E AK6 I 0.0 r DELTAE I DNSOIL * CPFUSE * AMOIST * XINCRE * YINCRE * +DELIl(I,L+1) / TIMSTP CHOICE I 22 END IF C C if node is diagonal to the node that borders the bottom C corner of tube C ELSE IF (J.EQ.4 .AND. K.EQ.TUBEHT-2) THEN IF (HICE2(I,L) .LE. (0.5*ZINCRE)) THEN IF (HICE2(I,L+1) .GT. (0.5*ZINCRE)) THEN C there is some freezing occuring AHMAJ - TIMSTP / (AMASOL . CPSOIL) AKl I 0.0 AHZ I 0.0 AK3 I 0.0 AK4 I 0.0 AKS I 0.0 AKG I 0.0 DELTAE I DNSOIL * CPFUSE * AMOIST * XINCRE * +DELIZ(I,L+1) * (2.0 * (HICE2(I,L) - 0.5 * ZINCRE) + +DELIZ(I,L+1)) / TIMSTP CHOICE I 23 ELSE 1259 C no freezing occurs AKMAJ I TIMSTP / (AMASOL * CPSOIL) AKl I AKSOLl I XINCRE * ZINCRE / YINCRE AKZ I AKSOLZ * XINCRE * ZINCRE / YINCRE AK3 I AKSOLJ * XINCRE * YINCRE / ZINCRE AK4 I AKSOL4 * XINCRE * YINCRE / ZINCRE AKS I 0.0 AK6 I 0.0 DELTAE I 0.0 CHOICE I 24 END IF C ELSE C there could be some freezing or defrosting AKMAJ I TIMSTP / (AMASOL * CPSOIL) AKl I 0.0 AK2 I 0.0 AK3 I 0.0 AK4 I 0.0 AKS I 0.0 A36 I 0.0 DELTAE I DNSOIL * CPFUSE * AMOIST * XINCRE * +DEL12(I,L+1) * (2.0 . (HICE2(I,L) - 0.5 . ZINCRE) + +DELIZ(I,L+l)) / TIMSTP CHOICE I 25 END IF if node is directly to the right of the node that borders the bottom corner of tube 0000 ELSE IF (J.EQ.4 .AND. X.EQ.TUBEHT-l) THEN IE (HICE2(I,L) .LE. (0.5*ZINCRE)) THEN IF (HICE2(I,L+1) .GT. (O.5*ZINCRE)) THEN there is some freezing occuring AKMAJ I TIMSTP / (AMASOL * CPSOIL) AKl I 0.0 0 2130 AKS I 0.0 AK6 I 0.0 DELTAE I DNSOIL * CPFUSE * AMOIST * XINCRE * ZINCRE * +DELI2(I,L+1) / TIMSTP CHOICE I 26 ELSE C no freezing occurs AKMAJ - TIMSTP / (AMASOL * CPSOIL) AKl I AKSOLl * XINCRE * ZINCRE / YINCRE AK2 I AKSOLz * XINCRE * ZINCRE / YINCRE AK3 I AKSOLJ * RINCRE * YINCRE / ZINCRE AK4 I AKSOL4 * XINCRE * YINCRE / ZINCRE AKS I 0.0 AKG I 0.0 DELTAE I 0.0 CHOICE I 27 END IF C ELSE C there could be some freezing or defrosting AKMAJ - TIMSTP / (AMASOL . CPSOIL) AKl I 0.0 AKZ I 0.0 AK3 I 0.0 AK4 I 0.0 AKS .‘ 0.0 AK6 I 0.0 DELTAE I DNSOIL * CPFUSE * AMOIST * XINCRE * ZINCRE * +DELI2(I,L&l) / TIMSTP CHOICE I 28 END IF 131 C if node is directly to the right of the node that borders C the center of side face of tube C ELSE IF (J.EQ.4 .AND. X.EQ.TUBEHT) THEN IF (HICEl(I,L) .LE. (0.5*ZINCRE)) THEN IE (HICE1(I,L+1) .GT. (O.5*ZINCRE)) THEN C there is some freezing occuring AKMAJ I TIMSTP / (AMASOL * CPSOIL) ARI I 0.0 AKZ I 0.0 AK3 I 0.0 AK4 I 0.0 AKS I 0.0 AK6 I 0.0 DELTAE I 0.5 * DNSOIL * CPFUSE * AMOIST * XINCRE * +SQRT(ZINCRE ** 2.0 + ABS(HICE1(I,L) - HICE2(I,L)) ** 2.0) * +(DELIl(I,L+l) + DELI2(I,L*1)) / TIMSTP CHOICE I 29 C ELSE C no freezing occurs AKMAJ I TIMSTP / (AMASOL * CPSOIL) AXl I AKSOLl * XINCRE * ZINCRE / YINCRE AK2 I AKSOLZ * XINCRE * ZINCRE / YINCRE AK3 I AKSOLJ I XINCRE * YINCRE / ZINCRE AK4 I AKSOLA * XINCRE * YINCRE / ZINCRE AKS I 0.0 AK6 I 0.0 DELTAE I 0.0 CHOICE I 30 END IF ELSE C there could be some freezing or defrosting AKMAJ - TIMSTP / (AMASOL w CPSOIL) AKl I 0.0 2132 AH2 - 0.0 AH3 - 0.0 AH4 - 0.0 AKS - 0.0 AK6 - 0.0 DELTAE - 0.5 * DNSOIL . CPFUSE * AMOIST 9 XINCRE * +SQRT(ZINCRE *9 2.0 + ABS(HICE1(I,L) - HICE2(I,L)) ** 2.0) * +(DELIl(I,L+l) + DELIz(I,L+1)) / TIMSTP CHOICE - 31 END IE if node is directly to the right of the node that borders the top corner of tube ' 0000 ELSE IF (J.EQ.4 .AND. K.EQ.TUBEHT+l) THEN IF (HICE1(I,L) .LE. (0.5*ZINCRE)) THEN IF (HICE1(I,L+1) .GT. (O.S*ZINCRE)) THEN C there is some freezing occuring AKMAJ I TIMSTP / (AMASOL * CPSOIL) AKl I 0.0 AKZ I 0.0 AK3 I 0.0 AK4 I 0.0 ANS I 0.0 AK6 I 0.0 DELTAE I DNSOIL * CPFUSE * AMOIST * XINCRE * ZINCRE * +DELIl(I,L+l) / TIMSTP CHOICE I 32 ELSE C no freezing occurs AKMAJ I TIMSTP / (AMASOL * CPSOIL) AKl I AKSOLl * XINCRE * ZINCRE / YINCRE AKZ I AKSOLZ * XINCRE * ZINCRE / YINCRE AKJ I AKSOL3 * XINCRE * YINCRE / ZINCRE AK4 I AKSOL4 * XINCRE * YINCRE / ZINCRE 133 AHS I 0.0 AK6 I 0.0 DELTAE I 0.0 CHOICE I 33 END IF ELSE C there could be some freezing or defrosting AKMAJ I TIMSTP / (AMASOL * CPSOIL) AKl I 0.0 AKZ I 0.0 AK3 I 0.0 AK4 I 0.0 AKS I 0.0 AKS I 0.0 DELTAE I DNSOIL * CPFUSE * AMOIST * XINCRE * ZINCRE * +DELIl(I,L+l) / TIMSTP CHOICE I 34 END IF C C if node is diagonal to the node that borders the top C corner of tube C ELSE IF (3.20.4 .AND. K.EQ.TUBEHT+2) THEN IF (HICEl(I,L) .LE. (O.5*ZINCRE)) THEN IF (HICEl(I,L+1) .GT. (O.5*ZINCRE)) THEN C there is some freezing occuring AKMAJ I TIMSTP / (AMASOL * CPSOIL) AKl I 0.0 AKZ I 0.0 AK: I 0.0 AK4 I 0.0 AKS I 0.0 AKG I 0.0 DELTAE I DNSOIL * CPFUSE * AMOIST * XINCRE * +DELIl(I,L+1) * (2.0 * (HICEI(I,L) - 0.5 * ZINCRE) + 134’ +DEL11(I,L+1)) / TIMSTP CHOICE I 35 C ELSE C no freezing occurs AKMAJ I TIMSTP / (AMASOL * CPSOIL) AX1 I AKSOL1 * XINCRE * ZINCRE / YINCRE AKZ I AKSOL2 * XINCRE * ZINCRE / YINCRE AK3 I AHSOL3 * XINCRE * YINCRE / ZINCRE AK4 I AKSOL4 * XINCRE * YINCRE / ZINCRE AKS I 0.0 AK6 I 0.0 DELTAE I 0.0 CHOICE I 36 END IF ELSE C there could be some freezing or defrosting AKMAJ I TIMSTP / (AMASOL * CPSOIL) ARI I 0.0 AKZ I 0.0 AK3 I 0.0 AK4 I 0.0 AKS I 0.0 AK6 I 0.0 DELTAE I DNSOIL * CPFUSE * AMOIST * XINCRE * +DELIl(I,L+1) * (2.0 I (HICE1(I,L) - 0.5 * ZINCRE) + +DELIl(I,L+l)) / TIMSTP CHOICE I 37 END IF C C if node borders the ground surface C ELSE IE (x .EQ. zNUM+1) THEN AKMAJ - TIMSTP / (AMASOL * CPSOIL) AKl - AHSOLI . XINCRE . ZINCRE / YINCRE 135 AKZ I AKSOL2 * XINCRE * ZINCRE / YINCRE 7 3 u I 2.0 * AKSOL3 * XINCRE * YINCRE / ZINCRE AK4 I 2.0 * HSUR * XINCRE * YINCRE AKS I 0.0 AK6 I 0.0 DELTAE I 0.0 CHOICE I 38 C if node is the center of the tube ELSE IF (J.EQ.2 .AND. K.EQ.TUBEHT) THEN C soil temperature need not be calculated here - proceed to AKMAJ I TIMSTP / (AMASOL * CPSOIL) C next element GO TO 340 C ELSE C C for other interior nodes C ARI I AKSOLl AK2 I AKSOLZ AK: I AKSOLJ AK4 I AKSOL4 ANS I 0.0 AKG I 0.0 DELTAE I 0.0 CHOICE I 39 END IF C i i i * XINCRE * ZINCRE / YINCRE XINCRE * ZINCRE / YINCRE XINCRE * YINCRE / ZINCRE XINCRE * YINCRE / ZINCRE ,C i****i**************************i**i********ii*&******************** C calculate the temperature of the element C ***********tittti****i********************************************** C TEMP(I,J,K,L+1) - AW . 2136 + (ARI 9 (TEMP(I,J-1,R,L) - TEMP(I,J,R,L)) + + AR2 9 (TEMP(I,J+I,R,L) - TEMP(I.J.K.L)) + + AR3 9 (TEMP(I,J,R-1,L) - TEMP(I,J,K,L)) + + AK4 9 (TEMP(I,J,R+I,L) - TEMP(I,J,R,L)) + .9 ARE 9 (TEMP(I,J-1,TUBEHT,L) - TEMP(I,J,R,L)) + + AKS 9 (TEMP(I,J-1,TUBEHT,L) ~ TEMP(I,J,R,L)) + + DELTAE) + TEMP(I,J,R,L) C C ***************i**************ii***********************************t C calculate ice formations and make temperature adjustments C it****************************************************************** C . IF (CHOICE .E0. 1) THEN IE (TEMP(I,J,R,L+1) .LT. 32.0) THEN IE (HICEZ(I,L) .LE. (O.5*ZINCRE)) THEN AKMAJZ - TIMSTP / (2.0 9 DNSOIL 9 CPFUSE 9 AMOIST 9 +xINCRE 9 YINCRE) AROO - AMASOL 9 CPSOLC / TIMSTP DELIz(I,L+1) - AKMAJZ 9 + (AROO 9 (TEMP(I,J,R,L+I) - TEMP(I,J,R,L))) IE (DELIz(I,L+1) .LT. 0.0) THEN IF (HICE2(I,L) .EQ. 0.0) THEN DELIZ(I,L+1) - ABS(DELIZ(I,L+1)) HICE2(I,L+l) . HICE2(I,L) + DELIz(I,L+1) ELSE HICE2(I,L+1) - HICE2(I,L) + DELIz(I,L+1) IE (HICE2(I,L+1) .LT. 0.0) THEN HICE2(I,L+1) . 0.0 END IE END IE ELSE HICE2(I,L+1) - HICE2(I,L) + DELIZ(I,L+l) END IE 137' IE (HICE2(I,L+1) .LT. (0.59zINCRE)) THEN TEMP(I,J,R,L+I) - 32.0 END IF END IF ELSE HICE2(I,L+1) I HICE2(I,L) END IF ELSE IF (CHOICE .EQ. 2) THEN IF (HICE2(I,L) .GT. (O.5*ZINCRE)) THEN AKMAJZ I TIMSTP / (DNSOIL * CPFUSE * AMOIST * XINCRE * +YINCRE) AKOO I AMASOL 9 CPSOIL / TIMSTP DELIZ(I,L+1) I AKMAJZ * + (AKOO * (TEMP(I,J,K,L+1) - TEMP(I,J,K,L))) HICE2(I,L+1) I HICE2(I,L) + DELI2(I,L+1) IE (HICE2(I,L+I) .GT. (O.5*ZINCRE)) THEN IF (HIC82(I,L+1) .LT. (I.592INCRE)) THEN IF (TEMP(I,J,R,L+I) .LT. 32.0) THEN TEMP(I,J,R,L+1) . 32.0 END IE ELSE WRITE(8,*)' ICE FORMATION Is GREATER THAN OR EQUAL + T0 3/2 ZINCRE’ GO TO 320 END IF END IF END IF 2138 ELSE IF (CHOICE .EQ. 3) THEN IF (TEMP(I,J,K,L+1) .LT. 32.0) THEN IE (HICEI(I,L) .LE. (0.59zINCRE)) THEN ARMAJz - TIMSTP / (2.0 9 DNSOIL 9 CPFUSE 9 AMOIST 9 +XINCRE 9 YINCRE) AROO - AMASOL 9 CPSOLC / TIMSTP DELIl(I,L*l) - ARMAJz 9 + (AROO 9 (TEMP(I,J,R,L+1) - TEMP(I,J,R,L))) IE (DELIl(I,L+l) .LT. 0.0) THEN IF (HICEI(I,L) .EQ. 0.0) THEN DELII(I,L+I) - ABS(DELIl(I,L+l)) , HICEI(I,L+1) - HICE1(I,L) + DELII(I,L+1) ELSE HICE1(I,L+I) - HICEl(I,L) + DELII(I,L+I) IE (HICEI(I,L+1) .LT. 0.0) THEN HICEI(I,L+1) - 0.0 END IF END IF ELSE HICEl(I,L+l) - HICE1(I,L) + DELIl(I,L+1) END IF IF (HICEI(I,Lwl) .LT. (0.5*ZINCRE)) THEN TEMP(I,J,R,L+I) a 32.0 END IF END IF ELSE HICE1(I,L+1) s HICE1(I,L) END IF ELSE IF (CHOICE .EQ. 4) THEN IF (HICEI(I,L) .GT. (0.5921NCRE)) THEN 2139 AKMAJZ - TIMSTP / (DNSOIL 9 CPFUSE 9 AMOIST 9 XINCRE 9 +YINCRE) AROO - AMASOL 9 CPSOIL / TIMSTP DELII(I,L+1) - ARMAJz 9 + (AROO 9 (TEMP(I,J,H,L+I) - TEMP(I,J,K,L))) C IE (HICEI(I,Lfil) .GT. (O.5*ZINCRE)) THEN IF (HICEI(I,L+1) .LT. (1.59zINCRE)) THEN IF (TEMP(I,J,R,L+I) .LT. 32.0) THEN TEMP(I,J,R,L+I) . 32.0 END IF ELSE WRITE(8,*)’ ICE FORMATION IS GREATER THAN OR EQUAL + T0 3/2 ZINCRE’ GO TO 320 END IF END IF END IF C C ELSE IE (CHOICE.EO.5 .OR. CHOICE.EQ.6 .OR. CHOICE.EQ.9 .OR. + CHOICE.EQ.11 .OR. CHOICE.EO.12 .OR. CHOICE.EQ.13 .OR. + CHOICE.EQ.16 .OR. CHOICE.EQ.17 .OR. CHOICE.EQ.20 .OR. + CHOICE.EQ.22) THEN C if temperature is less than 32 degrees, set temperature to 32 degrees IF (TEMP(I,J,K,L+1) .LT. 32.0) THEN TEMP(I.J.X,L+1) - 32.0 END IF C C . ELSE IF (CHOICE.EQ.23 .OR. CHOICE.EQ.25 .OR. CHOICE.EQ.26 .OR. + CHOICE.EQ.28 .OR. CHOICE.EQ.29 .OR. CHOICE.EQ.31 .OR. + CHOICE.EQ.32 .OR. CHOICE.EO.34 .OR. CHOICE.EQ.35 .OR. + CHOICE.EQ.37) THEN C if temperature is less than 32 degrees, set temperature to 32 degrees IE (TEMP(I.J.R,L91) .LT. 32.0) THEN 2140 TEMP(I,J,R,L+1) - 32.0 END IF END IF 340 CONTINUE 330 CONTINUE 320 CONTINUE ***********************fi******************************************** reset tube entrance temperature to current ambient temperature ******t***********ttttittiitiififiitt*tti***************************** 000000 TEMP(1,2,TUBEHT,L+1) I TATM ***********i******************************************************** considering the plane of symmetry around the tube in the y-z plane - fill temp(i,l,k,l+1) with temp(i,3,k,l+l) *************************************i****************************** 000000 DO 400 A I 2,XNUM+1 DO 410 B I 2,ZNUM+1 TEMP(A,1,E,L&1) I TEMP(A,3,B,L+1) 410 CONTINUE 400 CONTINUE C C i******************************************************************* C match end wall temperature to satisfy the insulation condition C ****************i*************************************************** C TEMP(1,1,TUEEHT,L+1) I TEMP(2,1,TUBEHT,L+1) 141. TEMP(1,3,TUBEHT,L&1) I TEMP(2,3,TUBEHT,L+1) TEMP(1,2,TUBEHTI1,L+1) I TEMP(2,2,TUBEHTI1,L+1) TEMP(1,2,TUBEHT+1,L+1) I TEMP(2,2,TUBEHT+1,L+1) C C ******i*iifiiiifi****fi***fi**i******it********************************* C calculate the temperature of air in the tube C tti*********************************tit**********************t****** C AKMAJl I XINCRE / (AMASAR * CPAIR * VELCTY) AH1 I 2.0 * H * XINCRE * ZINCRE AH2 I 2.0 * H * XINCRE * YINCRE DO 420 I I 2,XNUM+1 TEMP(I,2,TUBEHT,L+1) I AKMAJl * +(AH1 9 (TEMP(I-1,1,TUBEHT,L+1) - TEMP(I-l,2,TUBEHT,L+l)) + + AH1 * (TEMP(I-1,3,TUBEHT,L+1) I TEMP(II1,2,TUBEHT,L+1)) + + AH2 * (TEMP(I-l,2,TUBEHTI1,L+1) I TEMP(I-1,2,TUBEHT,L+1)) + + AH2 * (TEMP(I-l,2,TUBEHT+1,L+1) - TEMP(I-l,2,TUBEHT,L+l))) + + TEMP(II1,2,TUBEHT,L+1) 420 CONTINUE ***i************************************tit************************* test for temperature stability - if stable terminate program *************i**fi*iiififi*tti*****iti*i**i*iiiiiii******************** 00000 TDIFF I ABS(TEMP(XNUM+1,2,TUBEHT,L+1) I TEMP(XNUM+1,2,TUBEHT,L +)) IF (TDIEE .LT. TDIFFS) GO TO 800 ****i*******t****i***fi*§**fi*****************fi*********************** move 1+1 temperatures to l for next loop *t*************t**t************************************************* 00000 DO 430 A I 2,XNUM+1 DO 440 B I 1,YNUM+2 DO 450 C I 2,2NUM+1 TEMP(A,B,C,L) I TEMP(A,B,C,L+1) 142 450 CONTINUE 440 CONTINUE 430 CONTINUE C C *fiititiiii*iit**t*fl*t*fi*t*tiiii*iiiititiiititi*****i*************t** C move 1+1 ice formation to l for next loop 2 C tit***************************iiiititifiittittt********************** C DO 460 A I 2,XNUM+1 HICEl(A,L) I HICEl(A,L+l) HICE2(A,L) I HICEZ(A,L+1) 460 CONTINUE C C 310 CONTINUE C C ***********************i******************************************** C write out the temperature of the soil elements for every TIMLOP C ti*ttiiifiifiiiiiiiiiitiittttttfiiittt********************************* C WRITE(8,*)' TEMPERATURE OF NODES FOR EACH X-INCREMENT :' WRITE(8,10) DO 600 I I 2,XNUM+1 WRITE(8,50)I 50 FORMAT(1X,’ X I ’,I2) WRITE(8,*) WRITE(8,90) 90 FORMAT(1X,’K’,4X,’J I 2’,SX,’J I 3',5X,'J I 4',5X,’J = 5’, +5x,'J I 6',5X,’J I 7') DO 500 K I ZNUM+1,l,-l WRITE(8,40)K,TEMP(I,2,K,L),TEMP(I,3,K,L), +TEMP(I,4,K,L),TEMP(I,5,K,L),TEMP(I,6,K,L),TEMP(I,7,K,L) 40 FORMAT(1X,I2,6F11.7) 500 CONTINUE C wRITE(3,9) 143 WRITE(8,91) 91 FORMAT(1X,’K',4X,’J I 8’,5X,’J I 9',4X,’J I 10’,4X,’J = 11’ +,4X,’J I 12’,4X,’J I 13') DO 510 K I ZNUM+1,1,-1 WRITE(8,40)K,TEMP(I,8,K,L),TEMP(I,9,K,L), +TEMP(I,10,K,L),TENP(I,11,K,L),TEMP(I,12,K,L), +TEMP(I,13,K,L) 510 CONTINUE C WRITE(8,*) WRITE(8,92) 92 FORMAT(1X,'K',3X,'J - 14',4X,'J . lS’,4X,’J - 16',4X, +’J . 17') DO 520 K - ZNUM+1,1,-1 WRITE(8,41)K,TENP(I,14.x,L),TEMP(I,15,K,L), +TEMP(I,16,K,L),TEMP(I,17,K,L) 41 FORMAT(1X,IZ,4E11.7) 520 CONTINUE C WRITE(8,*) soo CONTINUE C c ******************************************************************** C write out the temperature of air in the tube for every TIMLOP C *****************ti*****t**i*****i******iii*iiii******************** C WRITE(8,*)’ TEMPERATURE OF AIR IN THE TUBE :' WRITE(8,10) DO 700 I I 2,11 WRITE(8,51)TEMP(I,2,TUBEHT,L+1),TEMP(I+10,2,TUBEHT,L+1), +TEMP(I+20,2,TUBEHT,L+1),TENP(I+30,2,TUBEHT,L+1), +TEMP(I+40,2,TUBEHT,L+1) 51 FORMAT(1X,5F11;7) 7oo CONTINUE C C ******************i**i*******fi****i*******fi************************* ‘144 C write out the error in temperature difference C **************t***tittttitttttttttiti*ttitiitttti**************tt*** _C WRITE(8,*) WRITE(8,52)TDIEE 52 FORMAT(1X,' ERROR IS ',F12.8) C ********************************i*********************************** C calculate time for every TIMLOP C *******************************************************************t C . HOURS I TIMSTP * TIMLOP * T WRITE(*,53)TEMP(XNUM+1,2,TUBEHT,L+1) WRITE(8,53)TEMP(XNUM+1,2,TUBEHT,L+1) 53 FORMAT(1X,’ OUTPUT TEMPERATURE OF TUBE IS ',F11.7) WRITE(*,54)HOURS WRITE(8,54)NOURS 54 FORMAT(1X,’ TIME IS ',F12.8,' HOURS’) C C ***i*********ifii********i**fi*****t********************************** C write ice formation for every TIMLOP C ***********************************t******************************** C WRITE(8,*) WRITE(8,*)’ THICKNESS OF ICE FORMED BELOW THE TUBE:’ DO 900 I I 2,11 WRITE(8,30)HICE2(I,L),HICE2(I+10,L),HICE2(I+20,L), +HICE2(I+30,L),HICE2(I+40,L) 30 FORMAT(1X,SF10.8) 900 CONTINUE C WRITE(8,*) WRITE(8,*)’ THICKNESS OE ICE FORMED ABOVE THE TUBE:’ DO 910 I - 2,11 WRITE(8,31)HICE1(I,L),HICE1(I+10,L),HICE1(I+20,L), +HICE1(I+30,L),HICE1(I+40,L) 31 910 C 'C 300 00000 00000 800 55 56 1145 FORMAT(1X,SF10.8) CONTINUE CONTINUE ****************************i******§******************************** if maximum time step is reached write out message ***ii************iiitttittittttfiiiiii*****i**fi********************** WRITE(8,*) WRITE(8,lO) WRITE(8,10) WRITE(*,*)' MAXIMUM TIME EXCEEDED' WRITE(8,*)’ MAXINUN TINE EXCEEDED' GO TO 810 **ii*iiiiifi******i*i***iiiiiitflfitiit*****i*************************i if stable temperature is achieved write out message ***t*****tt*itiiifiittiiifit*tiflttiiii******************************** WRITE(8,*) WRITE(8,10) WRITE(8,lO) WRITE(*,*)' OUTPUT TEMPERATURE IS STABLE' WRITE(8,*)’ OUTPUT TEMPERATURE Is STABLE' HOURS - (TIMSTP * (T-l) . TIMLOP) + (TIMSTP * P) WRITE(*,55)TEMP(XNUM+1,2,TUBEHT,L+l) WRITE(8,55)TEMP(XNUM+1,2,TUBEHT,L+1) FORMAT(1X,’ OUTPUT TEMPERATURE OF TUBE IS ',Fll.7) WRITE(*,56)HOURS WRITE(8,S6)HOURS FORMAT(1X,' TIME IS ',F12.8,’ HOURS’) WRITE(8,*)' TEHPERATURE OF AIR IN THE TUBE:' WRITE(8,lO) 750 760 770 810 820 146 DO 750 I - 2,11 ‘ WRITE(8,51)TEMP(I,2,TUBERT,L+l),TEMP(I+10,2,TUBEHT,L+1), +TEMP(I+20,2,TUBEHT,L&1),TEHP(I+30,2,TUBEHT,L+1), +TEMP(I+4O,2,TUEEMT,L+1) CONTINUE WRITE(8,*) WRITE(8,*)’ THICKNESS OE ICE FORMED BELOW THE TUBEz’ DO 760 I - 2,11 WRITE(8,30)HICE2(I,L#1),HICE2(I+10,L+1),HICE2(I+20,L+l), +HICE2(I+30,L+1),HICE2(I+40,L+1) CONTINUE WRITE(8,*) WRITE(8,*)’ THICKNESS OF ICE FORMED ABOVE THE TUBE:' DO 770 I I 2,11 WRITE(8,30)HICE1(I,L+1),HICE1(I+10,L+1),HICE1(I+20,L+1), +RICE1(I+30,L+1),HICE1(I+40,L+1) CONTINUE WRITE(8,*) GO TO 820 WRITE(*,*)’ CONTINUE? (YES . 1)' READ(*,*)ANS IF (ANS .EQ. 1) GO TO 20 WRITE(*,*)’ PROGRAM TERMINATED' WRITE(8,*)’ PROGRAM TERMINATED’ CLOSE(8) STOP END APPENDIX C DERIVATION OF HEAT TRANSFER EQUATIONS 147 Appendix C DERIVATION OF HEAT TRANSFER EQUATIONS There are several basic heat transfer equations which describe the heat transfer mechanism.in the models. One-Dimensional Model 1. Energy balance on an element of air within the tube, assuming steady state is as follows (derivation of equation 3.1): m C a a k _ k = k _ k _K§E- u (Ti+1,0 Ti,0) h A (Ti,1 Ti,0) or k k _ k _ i+1,o ' T1,o) ‘ h znrle (Ti,1 Ti,o) Solving for T§+1 0 : I k 2h Ax k k k i+1,o = pacparlu (Ti,1 - Ti,o) + T1,0 148 2. Energy balance on the soil element is as follows (derivation of equation 3.5): Ar m c 2n(r - -—)Ax 8 ps k+1 _ = j 2 k _ k At (Ti,j Tk,j) ksl Ar (Tl,j-1 Ti,j) 2n(r. - %£)Ax k k +1: 3 (T. . T. ) 52 Ar 1,]+1 1,j Solving for T¥+1 : 1.1 Ar Tk+i = —AE—— [ k 2n(rj §_)Ax (T3 . - T3 .) 1,3 mscps 31 Ar 1,3-1 1,3 2n(r. + 95M): + k 3 2 (T1? . Tk ) ] 82 Ar 1IJ+1 ll] k + . T1,: Here, heat transfer in the x-direction is taken to be negligible since it takes so long for temperature to Change with x. 3. When the node borders on the tube, the equation for the soil element is used, as in number 2, with the substitution of the approprite subscripts and application of the boundary condition (derivation of equation 3.2) . The equation used as in number 2 is as follows: 149 m c 2n(r - )Ax s s k+1 k l 2 k k XE (Ti,1 Ti,1) ks1 Kr (T1,O' T1,l) Ar + k hurl 2—)Ax (Tk - Tk ) (C 1) 32 Ar i,2 1,1 ' where T]; 0,is an image point. Applying and numerically translating the boundary condition given by equation 2.4 results the following: Tk - Tk 1,2 1,0' 3 k k ksA 2Ar hA (Ti,1 Ti,o) or k k T - T. _ Ar 1,1 1,0’ kslzn(r1 §—)Ax 2Ar + Tk - TR A; i,2 1,; = k _ k k822n(r1 + 2 )Ax 2At h 2nr1Ax (Ti,0 Ti,1) Solving for Tk , : 1,0 TI,0' = Ar Ar Eh 4nr1Ax (TIE,0 - T§,1) k812n(r1 - §—)Ax 2n(r1+ %£)Ax k k r ,2 1,1 + T‘.‘ 150 Substituting this image point into equation C.l, the following equation is thus obtained: m C s s k+1 k _ Tk _ Ar 2n(r1+ §—)Ax + 2ks2 Ar where m5 = psAvelement Ar 2 2 psn [(r1+ §—) - r1 ]Ax = psn (r1+ %£)ArAx k+1 i,l Solving for T Tlifii 3 Air [11 “"1“" (T); o ’ T]; 1) ' — I I psn(r1+ 4 )ArAxcps Ar + k 4n(r1+ §—)Ax (T3 - Tk ) s2 Ar 1,2 1,1 + Tli,1 4- When ice formed at the boundary of the tube is less than Ar/z, the derivation of equation 3.3 follows: 151 Ar psn(r1+ z—)ArAxcps (Tk+1 - Tk ) At 1,1 1,1 = Release of Latent Heat where, Ar. Release of Latent Heat = psLM 2nR Ax -Z%E§ Solving for Arice: c 3 At 3 25 A; k+1 _ k Arice psLM 2nR Ax[é t n(r1+ 4 )ArAx (Ti,l Ti,1)] 5. When ice is formed at the next soil element, the derivation Of equation 3.4 follows: mscps k+1 k At (Ti,2 - Ti,2) = Release of Latent Heat where, ms: psAvelement Ar 2 _ _ Ar 2 = 951! [(r2+ 2—-) (r2 5‘) 1 Ax = p82nr2ArAx Ar. ice Release of Latent Heat = psLM 2nR Ax At Solving for Arice: 152 k+1 Arice = E_Efi_%fi§—K§ [——;E— 2nr2 ArAx (Ti - Ti,2)] Two-Dimensional Model 1. The steady state energy balance on the air element in the tube is as follows (derivation of equation 3.6): 1 Solving for T1+1 O o: I I Ti+1,o,o a fi;%:;fi [ 2h + 2h + 2h + 2h + T1,0,0 AxAz AxAz AxAy AxAy AXAZ (Ti,-1,o ' T1,0,0 AXAZ (Ti,1,o ' Ti,o,o) AXAY (Ti,o,-1 ' Ti,o,o) AXAY (Ti,o,1 Ti,o,o) (Ti,-1,o ' Ti,o,o’ (Ti,1,o Ti,o,o) (Ti,o,-1 ' Ti,o,o) (Ti,o,1 Ti,o,o) ] 153 2. Energy balance on the soil element is as follows (derivation of equation 3.20): 22:55 (Tit%,k ‘ Ti,j,k) ‘ kslé§$5 (Ti,j-1,k ' Ti,j,k) + kszé%%£ (Ti,j+1,k ' Ti,j,k) + ksaé%%x (Ti,j,k-1 ' Ti,j,k) + ks4é%%x (Ti,j,k+1 ‘ Ti,j,k Solving for Ti+; k: . . Tit;,k ' figéfig [kslé%$£ (Ti,j-l,k ’ Ti,j,k) + kszé%%£ (Ti,j+1,k ' Ti,j,k) + kssé§gz (Ti,j,k-1 ' Ti,j,k) + ks4é%%x (Ti,j,k+1 ' Ti,j,k)] + Ti’j’k As in the one-dimensional model, heat transfer in the x-direction is also neglected since it takes so long for temperature to change with x. 3. When the node borders on the center of the bottom face of the tube, the equation obtained in number 2 is used. with the proper substitution of the appropriate subscripts 154 and boundary condition (derivation of equation 3.7) . The following equation as in number 2 is obtained: me 3 ps 1+1 _ l ' AxAz l _ 1 At (Ti,o,-1 Ti,o,-1) ksl-Ay— (T1,-1,-1 Ti,o,-1) AxAz l l + ksz‘Zy' (Ti,1,-1 Ti,o,-1) AxA 1 1 + ksz'KE (Ti,o,-2 Ti,o,-1) AxA 1 1 + karfix (Ti,o',o' Ti,o,-1) (C°2) where Ti 0’ 0' is an image point. Applying and numerically I I translating the boundary condition given by equation 2.8, the following equation results: T1 _ T1 _ 1,0',o' i,o,-2 = 1 _ 1 Rs 2Az h (Ti,O,-1 Ti,0,0) or 1 1 1 1 _k Ti,0',0' T1,o,-1 _ k T1,o,-1 Ti,o,-2 s4 2Az s3 2Az _ 1 1 h (Ti,o -1 T1,o,o) 1 . Solving for Ti,0’,0" 1 _ 2h A2 1 _ 1 T1,o',o' (Ti,o,o Ti,o,-1) s4 155 ks3 l l 1 + E34 (Ti,o,-2 Ti,O,-1) + Ti,O,-1 Substituting the image point into equation C.2, the following equation is thus obtained: m C 5 ps 1+1 _ l = AxAz l _ 1 At (Ti,o,-1 T1,o,-1) ksl-Ay— (Ti,-1,-1 Ti,o,-1) AxAz l 1 + ksz'Ay" (T1,1,-1 ’ Ti,o,-1) AxA 1 1 + 2ks3‘Z‘z‘x (Ti,o,-2 ' Ti,o,-1) 1 1 Solving for T1+1 : 1,0,-1 1+1 At AxAz l 1 T1,o,-1 a mscps [ks1 Ky (T1,-1,-1 Ti,o.-1) AxAz l l + ksz'Zy“ (Ti,1,-1 Ti,o,-1) AxA l l + sta'KEx (T1,o,-2 ' Ti,o,-1) 1 l + 2h AxAy (Ti,o,0 Ti,o,-l)] 1 + Ti,o,-1 4. When ice formed below the tube is less than Az/2, the derivation of equation 3.8 is as follows: 156 mscpc 1+1 1 At (Ti 0 _1 - Ti,O,-1) = Release of Latent Heat where, Ahic 2 Release of Latent Heat = ZpSLH AxAy __AEE— Solving for Ahicez: 1+1 1 Ahice2 2p LHt x y #[ Ktc(Ti1Ti,O,-l)]- 5. Energy balance on the soil element at the bottom corner of the tube (derivation of equation 3.10) is as follows: Internal 3 Conduction Heat Total Convection Distributed ,Rates along Four Rate along Two Source of Different Lanes Half-Lanes in Energy from Neighboring the y and z Nodes directions or 3m c s 5 1+1 _ l = AxAz _ 1 4At (Ti,1,-1 Ti,1,-1) ks1‘23" (Ti, 0, -1 Ti,1,-1) AxAz 1 + k32"— Ay (T1, 2, -1 ' T1,1,-1) AxA l + ks3 Az (Ti,1,-2 - Ti,1,-1) 157 AxA l 1 + R84 2 z (Ti,1,o ' Ti,1,-1) Az l l + h Ax -2--(Ti’o'o Ti,1,-1) Ay 1 _ 1 Solving for T1 ° 1,1,-1' 1+1 2At AxAz l TL1,-1 SE;E;; [ksl 3“" (Ti, 0, -1 Ti,1,-1) AxAz l + 2k52-2!" (TL 2,-1 ' Ti,1,-1) + 2ks3_AEX (Ti,1,-2 ‘ Ti,1,-1’ + ks4AEAX (Ti,1,o ’ Ti, 1, -1) + h AxAz (Ti 0 0 - Ti'1'_1) + h AxAy (Ti'o’0 - Ti’1'_1)] + Ti,1,-1 6. Energy balance on the soil element at the bottom corner Of the tube when ice formed is less than Az/2 (derivation Of equation 3.11) is as follows: 3m c s 3 1+1 1 4 t (Ti,l _1 Ti,l,-1) Release of Latent Heat where, Release of = p LH Ax (Ax + _5 + 2h + ) Ahice2 Latent Heat 2 2 ice2 hice2 Kt 158 . 1+1 , Solv1ng for Ti,1,-1' 1+1 = 2At Ti,1,-1 3m c [ psLM Ax (Ay + AZ 8 ps + 4h Ah ice2 ice2 + 2Ahice2) At ] "71111111111111?