““43“ LIBRARY
Michigan State
University l

 

This is to certify that the

thesis entitled
The Thermal Performance

Of An Underground Air Pipe
presented by

Brian Michael Leary

has been accepted towards fulfillment
of the requirements for
Master's Mechanical

degree in ——Eng—i—neering

Major professor

 

Date 4-24-86

 

0-7639 MS U is an Affirmative Action/Equal Opportunity Institution

i.

 

 

 

RETURNING MATERIALS:
)V1531u1 Place in book dron to
LJBRARJES remove this checkout from
.AuuumynaauL your record. FINES will

 

 

be charged if BEBE-is
returned after the date
stamped below.

 

 

 

 

THE THERMAL PERFORMANCE
OF AN UNDERGROUND AIR PIPE

By

Brian Michael Leary

A THESIS

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

MASTER OF SCIENCE
Department of Mechanical Engineering

1984

ti

pr
of

tc
to

to

of
to

of

 

ABSTRACT

THE THERMAL PERFORMANCE
OF AN UNDERGROUND AIR PIPE

By

Brian Michael Leary

The heat output of a typical air-to-air heat pump decreases as
the ambient temperature decreases. A potential solution to this
problem is to preheat the air supplied to the outside heat exchanger
of the heat pump by drawing it through an underground duct.

A finite difference analysis of the heat transfer from the soil
to an underground pipe in which air is flowing is performed in order
to determine the feasibility of this technique. The latent heat due
to the formation of ice lenses and the migration of moisture toward
the freezing front are included in the analysis. The results show
that this would be a viable technique for preheating the air supply
of an air-to-air heat pump significantly improving its capability
to supply heat at low ambient temperature while raising the coefficient

of performance.

.uuv

R IV:

;-‘

lxl.

rbe

.hu

5 l

ACKNOWLEDGEMENTS

I would like to thank Dr. Merle C. Potter for his help and
guidance throughout the preparation of this work. Thanks also
to Dr. J. Beck and C. R. St. Clair, Jr. for their helpful comments
and suggestions.

Additionally, I would like to thank Ms. Judy Duncan for her
patient and careful preparation of the final manuscript, Mrs. Deann
Heath for preparing many of the figures, and Mr. Paul Zang for his
help with generation of the figures detailing the results.

Finally, I would like to thank my parents, John and Patricia
Leary, and Darlene Bloomquist for their support and encouragement

throughout my academic career.

ii

 

 

f‘l

TABLE OF CONTENTS

Page
LIST OF TABLES ................................................. v
LIST OF FIGURES ................................................ vi
NOMENCLATURE ................................................... vii
CHAPTER 1 -- INTRODUCTION ...................................... l
1.1 Modifications to Improve Heat Pump
Performance .......................................... 3
1.2 Surface and Well Water Systems ....................... 5
1.3 Recirculating Systems ................................ 6
1.4 Recirculating Air Systems ............................ 9
1.5 Potential for Thermal Reservoir Systems .............. 10
1.6 The Buried Duct System ............................... 11
CHAPTER 2,-- ANALYSIS .......................................... 17
2.1 Determination of Energy Stored in
Soil Column .......................................... l7
2.2 Distance to which Heat Extraction
Affects the Surrounding Soil ......................... 18
2.3 Governing Heat Transfer Equation ..................... 19
2.4 Modes of Heat Transfer ............................... 20
2.5 Boundary Conditions .................................. 21
2.6 Initial Soil Temperature Distribution ................ 24
2.7 Soil Thermal Properties .............................. 24
2.8 Soil Conductivity and Heat Capacity .................. 24
2.9 Moisture Attraction .................................. 26
2.10 Nature of Air Flow in the Duct ....................... 26
2.11 Numerical Stability .................................. 30
CHAPTER 3 -- NUMERICAL MODEL ................................... 31
3.1 Soil Properties and Initial Temperatures ............. 31
3.2 Soil Temperature Profile ............................. 34
3.3 Determination of Time Step ........................... 34
3.4 Moisture Migration Function .......................... 36
3.5 Determination of Soil Heat Capacity
and Conductivity as Functions of
Temperature .......................................... 36
3.6 Derivation of Heat Transfer Equations ................ 37
3.7 Numerical Heat Transfer Equation ..................... 41
3.8 Soil State Change Adjustment ......................... 45
3.9 Preparation for Subsequent Time Step ................. 46

TABLE OF CONTENTS (continued)

Page
3.10 Adjustment for Insulation Condition
and Plane of Symmetry ............................... 45
3.11 Reinitialization of Soil Temperature Array .......... 50
3.12 Programming Technique ............................... 50
3.l3 Computers Used ...................................... 51
CHAPTER 4 -- SUMMARY OF RESULTS ............................... 52
4.1 System Feasibility .................................. 52
4.2 Economic Feasibility ................................ 50
CHAPTER 5 -- CONCLUSION ....................................... 62
CHAPTER 6 -- RECOMMENDATIONS .................................. 63
APPENDIX A -- HEAT TRANSFER PROGRAM ........................... 67
APPENDIX B -- INTERACTIVE MICROCOMPUTER PROGRAM ............... 78
APPENDIX C -- ONE PASS TEST PROGRAM ........................... 87

BIBLIOGRAPHY .................................................. 103

LIST OF TABLES

Table Page
1 Initial Values ............................................ 33
2 Initial Soil Temperature Distribution ..................... 35

3 State Change Adjustment Table ............................. 47

LIST OF FIGURES

Figure Page
1 Balance Point of a Heat Pump ............................ 2
2 Heat Only Heat Pump Performance ......................... 4
3 Well Hater Reinjection System ........................... 7
4 Buried Refrigeration Tubing Loop ........................ 8
5 Vertical Recirculating System ........................... 10
6 Air Recirculation in Crawl Space ........................ I]
7 Crawl Space One Pass and Bypass Configurations .......... I3
8 Proposed Air Pre-conditioning System .................... I4
9 Plane of Symmetry Bisecting Pipe Vertically ............. 22

10 Insulated Y-Z Planes at Ends of Pipe .................... 23
11 Soil Temperature Profile ................................ 25
12 Soil Thermal Capacity as a Function

of Temperature .......................................... 27
13 Ice Formation Around Buried Duct ........................ 28
14 Velocity Profile in Duct ................................ 29
15 Numerical Scheme ........................................ 32
16 Heat Balance on Soil Elements Surrounded

by Soil Elements ........................................ 38
17 Heat Transfer to 3 Soil Element which Borders

a Fluid Element ......................................... 40
18 Heat Transfer to Fluid Elements ......................... 42
19 Satisfaction of Insulated End Conditions ................ 48
20 Plane of Symmetry ....................................... 49
21 Outlet Temperature Versus Time at Four

Velocities .............................................. 53
22 Temperature Versus X-Position for Various

Times when the Air Velocity is One Foot per Second ...... 54
23 Temperature Versus X-Position for Various

Times when the Air Velocity is Two Feet per Second ...... 55
24 Temperature Versus X-Position for Various

Times when the Air Velocity is Four Feet per Second ..... 55
25 Temperature Versus Time at Various X-Positions in

the Tube when the Air Velocity is One Foot per

Second .................................................. 57
26 Temperature Versus Time at Various X-Positions in

the Tube when the Air Velocity is Two Feet per

Second .................................................. 58
27 Temperature Versus Time at Various X-Positions in

the Tube when the Air Velocity is Four Feet per

Second .................................................. 59
28 Potential Duct Configurations .................. - ......... 55

vi

NOMENCLATURE

area, FT2

thermal diffusivity, FT2 /s

heat capacity, BTU/1b- F

energy, FT-lb 2
convective heat transfer coefficient, BTU/FT
distance FT 2

thermal conductivity, BTU- FT/HR- FT -°F
distance, FT

moisture content, %

mass, 1b 3

density, 1b/FT

heat, BTU

e Reynolds number

temperature, °F

time, HR 2
overall heat transfer coefficient, BTU/FT -°F
velocity

kinematic viscosity, FT2 /5

change

differential

Zr-xxa'mnw D

artesian coordinates
Cartesian X
Cartesian Y
Cartesian Z

rray indices

Cartesian Z direction
Cartesian X direction
Cartesian Y direction
time

ubscripts

air

deep soil region
final

initial

soil

time

zom—h-hCLmM PRC—4H):- N-<><O QD<<C¢+HJUDUE

-°F

CHAPTER 1
INTRODUCTION

In recent years a large number of air-to-air heat pumps have been
installed in northern climates, primarily to provide the necessary
cooling but also to contribute to the heating load. Air-to-air heat
pumps are up to three times more energy efficient than electric re-
sistance heating (1). Additionally, because air-to—air heat pumps are
driven by electricity, they are less affected by shortages of, or
increases in, the prices of natural gas and oil. Unfortunately, even
though air-tO-air heat pumps save energy during periods of cool weather,
their use in northern climates, particularly in cold weather, leads to
a number of problems.

In extremely cold weather, the heat output Of an air-to-air heat
pump and the heat pump's coefficient of performance are drastically re-
duced. Eventually, as the temperature drops, the balance point of the
heat pump is reached, as depicted in Figure l. The balance point,
generally between 18 and 28 degrees Fahrenheit, is the point where the
increasing heat demand of the building equals the decreasing heat out-
put of the heat pump (2). Below this point supplemental heating is
needed, usually resistance heat. At ambient temperatures around 17
degrees Fahrenheit the typical air-to-air heat pump's coefficient of
performance is reduced to one and the maximum heat output has fallen so

low that electric resistance heating or a conventional furnace must be

ass; “mm: m eo ucwom mocmme .P umanm

Auwwccmcsau mwmcmoov menuecmaem»

mm 0m. 0

.ZEQQ mgcxixm

BDQHDQ QEDQ .DQI

35$
RZMEMQSOMQ hflwl .ZQI

 

used to supply the building's heat load. Even heat only heat pumps,
such as Janitrol's Nattsaver, designed specifically for cold weather,
suffer a large performance degradation. The Nattsaver's coefficient
Of performance relative to temperature is shown in Figure 2 (3).

Another problem which occurs if the air-to-air heat pump is used
to supply heat in severe weather is that the high demand on the com-
pressor causes frequent breakdowns (4). A problem of great concern to
utility companies is that since the air-to-air heat pump Operates two
to three times more efficiently than electric resistance heating at
moderate temperatures, when electric resistance heating is used to
supply the additional demand during severe weathen,a disproportionally
high peak load occurs. This forces utilities to maintain a large

portion Of their generating capacity in expensive peaking plants (5).

1.1 Modifications to Improve Heat Pump Performance

A recent Oak Ridge National Laboratory study (6), commissioned
specifically to determine areas of potential heat pump performance
improvement, concluded that an increase in the size of heat exchangers
and improvements in heat exchanger configuration, more efficient fans,
and improvements in compressor efficiency would produce important
increases in heat pump performance when the heat pump is operated in
moderate climates. However, none of these improvements significantly
affect the problems caused by heat pump operation at low temperatures.

One largely unexamined solution to the problem of improving the
performance of the air-to-air heat pump at temperatures below 30 degrees
is the possibility Of providing higher temperature air to the outside
heat exchanger. If a constant temperature thermal reservoir at 30

degrees Fahrenheit could be provided for the heat pump,the result

mocmscowcma gaze «no: z—co one: .N mmszu

.m ~35 VQMQSMH k ZMQQE v.
em I Q I Q 9 ON on 9v On 00

 

0.0

Q;

 

 

/ ON

/ so.

 

 

 

 

 

 

 

 

 

 

Om

 

would be higher energy efficiency, good compressor reliability and a
large reduction in peak electrical demands. An analysis based on the
distribution of heating degree days in Detroit, indicates that if the
effective operating range of the air-to-air heat pump could be extended
from an ambient temperature Of 35 degrees Fahrenheit to an ambient
temperature of 20 degrees Fahrenheit, the heat pump's contribution to
the heating load would be increased 50 percent; and if the heat pump
could operate efficiently down to an ambient temperature of 10 degrees
Fahrenheit, its contribution to the heating load would be increased 100

percent.

1.2 Surface and Well Water Systems

A number of methods Of providing a higher thermal reservoir
temperature to various types Of heat pumps have been tested. Surface
water has been used with water-tO-air heat pumps; however, surface
water availability, heat exchanger scaling, and surface water freezing
are serious problems in these systems when used in northern climates.
One of the first commercially successful approaches to providing
water-to-air heat pumps with a high quality thermal source in a
northern climate was to pump ground water from wells and pass the water
over the heat pump's outside heat exchanger. A system installed at
the Batelle Institute in Columbus, Ohio is supplied with 40 million
gallons of 54 degree Fahrenheit water each month from five, 16 inch
wells. The system has provided year-round coefficients Of performance
as high as 5.4 (7). The system at the Batelle Institute draws its
water from a sand and gravel aquifer that is fed by the Olentangy River

and discharges the water into the River after use. In other areas, the

problems of aquifer depletion and waste water disposal frequently
restrict or prevent the use of well water systems. In order to
eliminate these problems, water reinjection systems have been developed
which use a second well to return the water to the aquifer (Figure 3).
However, the high capital cost entailed in drilling a second well and
legal restrictions due to the possibility of ground water contamination

severely limit the use of reinjection systems (8).

1.3 Recirculating Systems

A number of systems have been designed in which a closed system
circulating a working fluid draws heat from soil. One of the first
such systems was designed and tested by British engineer John Sumner
(9). In his system refrigerant was circulated through 1250 feet Of
copper tubing buried in a 30 by 45 foot rectangle of soil (Figure 4).
The tubing constituted the external heat exchanger of the heat pump.
The heat pump then supplied heated water which was circulated beneath
the floor of the dwelling. The system successfully supplied the 84
million British Thermal Unit heating demand of the residence which it
served. The heat pump used was designed and built by Sumner and is
one of the only "soil-to-water" heat pumps in existance. (The re-
frigerant is considered the working fluid.) The high cost Of re-
frigeration quality tubing and the potential for expensive, environ-
mentally damaging refrigerant leaks precludes the large scale use of
such systems.

A similar system (10), tested at Brookhaven National Laboratory,
allowed a heat pump to supply the heating demand of an 1120 square
foot residence without supplemental heating at ambient temperatures

as low as minus 11 degrees Fahrenheit. The liquid-tO-aid heat

 

UMP [HEAT \ SURFACE
AEXCHANGER y

 

 

 

 

 

 

 

F _H

 

 

 

 

X\

 

 

 

SUPPLY PE/NJECT/ON
WELL WELL
M'\ /1V’\

 

 

AQU/FER

 

 

 

 

FIGURE 3. Well Water Reinjection System

P— 30 4.

m mREFR/GERATION
(JO/Ls

 

 

 

45

 

 

 

 

HEAT

PUMP ' 4

-L- [I

[— RESIDENCE 1

FIGURE 4. Buried Refrigeration Tubing

 

 

 

 

 

 

pump's external heat exchanger extracted the heat collected by a 507 foot
long coil of 1.5 inch medium density polyethylene pipe buried 4 feet deep
in which an antifreeze solution was circulated. The system's overall
seasonal coefficient of performance was 2.2. Unfortunately, horizontal
systems are virtually irreparable during the heating season due to the
difficulty inherent in excavating a trench several hundred feet long

in frozen ground.

Nater-to-air heat pumps have been coupled to vertical recirculat-
ing heat exchangers with great success. Several hundred systems
utilizing u-shaped tubes circulating brine, installed in 150 foot
deep wells, are in commercial and residential use in the Stillwater,
Oklahoma area (Figure 5). A University of Oklahoma study claims that
these systems are economically competitive with all conventional
sources of energy except price controlled natural gas (11). A sig-
nificant problem with all of the systems using recirculating liquids
is that leaks disable the system, and since they all require use of
an antifreeze solution, leaks could conceivably severely pollute

local ground water.

1.4 Recirculating Air System

A novel system using air recirculated through a crawl space was
recently tested at the University of Tennessee (Figure 6) (12). The
crawl space was thoroughly isolated from the rest Of the house by high
resistance insulation. The study revealed that conditioning the air
supplying the external heat exchanger of an air-to-air heat pump by
recirculating it through the crawl space beneath the residence provided

a 26.3 percent reduction in peak heating demand and a 7.7 percent

IO

\

RESIDENCE CHEAT PUMP
——-\

 

 

 

 

U—TUBE ‘—“‘

 

/50 '

 

 

 

FIGURE 5._ Vertical Recirculating System

 

 

—a
—a
HEAT PARTITION
pump

{——

9 REC/RCULAT/NG
AIR

 

 

 

 

 

 

 

a. Top View Crawl Space

 

 

DWELLING
INSULATION
I In; I
b. Side View

FIGURE 6. Recirculating Crawl Space Heat Pump System

 

 

12

reduction in peak cooling demand. Overall heating load was reduced
14.4 percent, while overall cooling load was increased 2.8 percent.
The reason for the increase in cooling load was that the system was
always Operated in the recirculation mode even when the crawl space

temperature was far in excess of the ambient temperature.

1.5 Potential for Thermal Reservoir Systems

The success of the systems discussed above illustrates the
tremendous potential that ground coupling a heat pump can have on its
performance. Unfortunately, in addition to the problems detailed above,
systems utilizing liquids cannot take advantage Of periods when heating
is required and the ambient temperature is higher than the temperature
of the earth source. In addition there is no simple, economical way
of retrofitting these systems to the many air-tO-air heat pumps already
in operation. The crawl space system has tremendous potential and its
performance could probably be significantly improved by designing it
to Operate in single pass mode and bypass modes when appropriate
(Figure 7). However, crawl spaces are not universally incorporated

into houses so this system does not have universal potential.

1.6 The Buried Duct System

A significant need exists for a system that can extract energy
from the ground, is universally available and that can be retrofitted
to existing air-to-air heat pumps. A system Of buried ducts through
which air is passed and then supplied to the heat pump satisfies all
the requirements (Figure 8). The system could be easily retrofitted

to existing air-tO-air heat pumps. It could be operated as a

 —1

13

 

J _—‘) AIR FLOW

 

 

HEAT ’
* PUMP

 

 

 

T )D)

 

 

a. One Pass Operation

 

E ill

HEAT

 

 

'_‘ PUMP

Til

 

 

 

 

 

 

b. Bypass Operation

FIGURE 7.

Crawl Space One PaSs and Bypass Operation

 

 

l4

Emumzm mcwcowuvucoonmca cw< ammoaoca .w mmaumd

 

 

 

 

 

 

 

 

 

 

 

 

T .3 L.
r... a.
.m
i _
$59me / 2‘ S \x/\J
$31 A E:

 

L\

\So e?

 

recir
smal'
tempt
DUO?
be d

heat

air
rat-
be

ten
eff
hee
tr:

thl

l5

recirculating system, as a one pass system or bypassed entirely. A
small system could be designed to moderate the affects of daily
temperature variations, a slightly larger system could sustain the heat
pump through periods of moderately cold weather; or a large system could
be designed to supply preheated air to the heat pump through the entire
heating system.

A preliminary numerical study (13) indicated that 0 degree ambient
air could be drawn through a 1 foot square, 50 foot long tube at a
rate of 1.39 feet per second and 38 degree Fahrenheit output air could
be supplied to the heat pump. The study did not include the convective
temperature terms, the heat contained in the mdisture in the soil, the
effect of moisture migration toward the subcooled soil, and the latent
heat of fusion released when soil moisture freezes, all of which in-
crease the amount of heat available to heat the air. Additionally,
the time period analyzed was less than one hour.

The potential of the buried duct system revealed by this study
merits additional detailed study. The Objective of this study Of the
buried duct system was to accurately estimate the heat transfer attain-
able over a significant period of time in the tube mentioned above.

To achieve this Objective a numerical simulation was developed. The
analysis simulates the heat transfer to a square duct from the surround-
ing soil using a finite difference method. Forward differences are

used in time and central differences are used in space. Parameters
which may be varied include soil type, soil thermal properties, soil
moisture content, air velocity in the duct, ambient temperature, and
deep soil temperature. The pipe size and length, tube depth beneath

the surface, and soil element size are also varied. The effect Of

l6

soil moisture content, latent heat of fusion of moisture in the soil,
the increase in thermal diffusivity of frozen soil, the increase in the
moisture content of cooling soil are included in the simulation. The
simulation was translated into the Fortran Five programming language
and a number of analyses were run on the CDC Cyber 750 and the Prime

7500 computers.

CHAPTER. 2
ANALYSIS

Determining the heat transfer from the soil to the air flowing
in a buried duct (Figure 8) involves a number of factors, some of which
may be difficult to account for due to the lack of experimental data.
Among the factors considered are the nature of the air flow in the
pipe, the convective heat transfer from the soil to the air, the heat
conduction through the soil, the effect of soil moisture on heat
capacity, moisture migration toward the duct, the freezing Of the soil
surrounding the upstream part of the duct, and changes in soil pro-

perties with temperature.

2.1 Determination of Energy Stored in Soil Column

It is important to determine the size of the soil block necessary
to supply the heat required for a building over an entire heating season.
The initial pipe burial depth was chosen as 5 feet and empirical re-
sults from liquid systems indicate that no heat transfer takes place
beyond 7.5 feet in the course of a heating season (14). Therefore, the
depth of the soil block from which heat would be extracted was selected
as 12.5 feet. A computer program was written which calculated the
amount of heat that could be extracted from a 1 foot square column of
soil extending 12.5 feet into the ground if it was cooled from its

initial temperature to a final temperature of 30 degrees Fahrenheit.

17

18

The vertical temperature distribution in the soil column was determined
by using the equation derived for the main program which assigned
temperatures between the average ambient temperature and the constant
deep earth temperature using a second order polynomial. The heat
capacity of the soil was calculated using the heat capacity Of dry
soil, plus the heat capacity Of the soil moisture, plus the latent
heat of fusion released during the freezing of the soil moisture. It
was found that this column Of soil could supply 40,232 BTU's of heat
and therefore a surface area of 2064 square feet would be required to
supply the 83,000,000 BTU heating requirement for a typical residence
(15). While the analysis assumed that the entire soil block froze,

it ignored heat transfer from outside the boundaries of the soil block

and therefore should be a conservative estimate.

2.2 Distance to which Heat Extraction Affects the Surrounding Soil
Another value which is required for the planning of a numerical
solution is the lateral distance into the soil to which the effect
of heat extraction would extend over the heating season. This can
be determined from both empirical results and dimensional analysis

using the formula

OL*t/L2 = .1
where
a = soil thermal diffusivity in ftZ/hr
t = time in hours
L = distance in feet.

19

Solving for L, where a = .00144 ft2/hr, and t = 3840 hours

1/2

P
I

- (at/.1)

(.00144*384O/.1)1/2 = 7.44 feet.

'—
II

This agrees with the results of an experimental study (16) performed
at the University of Oklahoma which indicated that the effects of heat
extraction would not extend further than 7.5 feet over an entire heating

season .

2.3 Governing Heat Transfer Equation
The governing heat transfer equation for this problem is the three

dimensional transient diffusion equation

2 2 2

3%*:—yi*§‘§=9§ii (2.1)
where

T = temperature in °F

c = specific heat in BTU/1b-°F

k = conductivity in BTU/FT-hr-°F

p = density in lb/FT3

In order to solve this equation, six boundary conditions and one

initial condition must be known. These are easily determined, since
the x, y, and 2 planes where the temperature can be assumed constant
are readily calculated and a reasonably accurate estimate of the

temperature distribution in the soil at the start of Operation of the
pipe can be made. However, the presence Of the pipe complicates the
situation. Neither the temperature nor the heat flux at the wall of

the pipe are known making it impossible to solve the equation

 

 

 

ar

0L

Wn

20

analytically. Therefore, it was necessary to solve this problem
numerically.

Using forward differences in time and central differences in
space, a numerical equation was produced which yielded results with an
error on the order Of the size Of the space step squared and the time

step.

2.4 Modes of Heat Transfer

Heat transfer between nodes takes place either through conduction
or convection. At the boundary between the top node and the atmosphere
and the boundary between the pipe nodes and the flowing air the con-

vection equation
q = UA(Ti-To) (2.2)

is used where h is the overall heat transfer coefficient from the soil

block to the air in the pipe (17).
Numerically this translates to:

q = UAXAY(Ti-To) (2,3)
for heat transfer in the z direction.

Everywhere else the heat transfer occurs as conduction through the soil;

the conduction equation for heat transfer in the z-direction used is
q = kAdT/dz (2.4)

where k represents the soil conductivity (18).

21

This is expressed numerically, for heat transfer in the z-direction,

as
q = (kAXAY/AZ)(AT)

where the derivative dT/dz has been approximated by AT/AZ and AT is
the temperature difference between neighboring soil elements.
Some heat is carried by moisture migrating toward the cooling

soil but this heat is small and is neglected.

2.5 Boundary Conditions

It is important to identify the area through which significant
heat transfer is taking place (Figure 9). Horizontally, the heat
extracted by the pipe should not cause significant temperature changes
beyond 7.5 feet from the tube since the part of the heating season
during which the temperature is below 30 degrees is less than 160 days
long and, as demonstrated in the previous section on the impact of
heat extraction, the effects will not extend beyond 7.5 feet in this
time period. Similarly, significant temperature changes would not be
present more than 7.5 feet below the pipe.

The heat transfer effects are symmetrical about a vertical plane
bisecting the pipe longitudinally. Therefore, the heat transfer need
be calculated on only one side Of the pipe. There would also be some
heat transfer into the soil block at each end Of the pipe, since the
pipe must pass vertically through the soil in these regions; modeling
this heat transfer is assumed small and is thus neglected. Therefore,

these planes were modeled as being insulated (Figure 10). This is a

22

mcozgcou @2258 use x322? to 28E

MES «ISSUE Sam ammoue J

.m was:

 

b——--———-——-_I---

 

.82; L

 

5233.209 BL

 

IE E 38 A

 

 

-d--——------q

tease n \
amnesia, to midi A A

23

xcecczom :o_um_=sem um mosey; Ni> uwum_:mc_

-r-

eQm

.op mmeHu

 

 

 

 

[/11 1 I

 

7 EEXHX Q

X

 

 

 

24

conservative assumption since there would be heat transfer from these

blocks to the region containing the horizontal pipe.

2.6 Initial Soil Temperature Distribution

The initial soil temperature distribution is described by a
second order polynomial with a sinusoid superimposed (Figure 11) (19).
The sinusoidal variation in temperature does not extend to a significant
depth since heat transfer occurs into the low temperature areas from
above and below as the fluctuations penetrate deeper into the soil,

tending to quickly eliminate the fluctuations.

2.7 Soil Thermal Properties

After the governing heat transfer equation is identified, the
initial and boundary conditions determined, and the geometry of
interest specified, the thermal properties Of the heat transfer medium
are considered. The thermal properties of soil depend upon composition,
partical size, and moisture content. Soil thermal properties are
empirically determined due to the enormous variation possible even in

similar soils.

2.8 Soil Conductivity and Heat Capacity

Some of the changes in soil thermal properties as functions Of
temperature are universal and can be specified. Soil, at temperatures
above 32 degrees Fahrenheit has relatively low thermal conductivity,
and relatively high heat capacity. Soil which is in the process of
freezing has a slightly better conductivity, and a very high apparent
heat capacity due to the release Of the latent heat of fusion of the

soil moisture as it freezes. The bulk of the water in the soil freezes

25

N

mrweogm mgzgmgmasoh Fvom .PF mmonm

FREQ

 

Mbfikmbm I N1

 

wtbkv‘mquxmk $2.3m“)? RZMQQDQT

 

mtbkv‘QMQIE 45% QMMQIIJ

 

EL

HMEQmm<kbtm

 

betwee
frozen
therna

capaci

2.9 M(

I
Of mo'
far as
funct'
increg
water
The mi
drOps
aSag
V3ile:

and v

26

between 30 and 32 degrees Fahrenheit (Figure 12) (20). Soil which is
frozen exhibits the highest thermal conductivity due to the high
thermal conductivity of ice relative to water, but a lower heat

capacity due to the lower heat capacity Of ice.

2.9 Moisture Attraction

A final important thermal effect present in soil is the phenomenon
of moisture attraction. As soil is cooled it draws moisture from as
far as six feet away (Figure 13) (21). The moisture increase is a
function of both time and temperature and both the rate of moisture
increase and the final moisture content is greatly enhanced if the
water table is less than three feet below the freezing front (22).
The moisture content in some soils doubles as the soil temperature
drops to freezing. The soil moisture content increase can be modeled
as a second order function Of temperature; however, the exact form
varies since the equations suggested are largely empirically derived

and vary greatly between researchers in the field (23).

2.10 Nature of the Air Flow in the Duct

The nature of the air flow in the duct is important to the heat
transfer. At a speed of 1 foot per second in a 12 inch duct the
Reynolds number of the flow (Re = 6250) is high enough to ensure com-
pletely turbulent flow, even near the duct entrance since large dis-
turbances occur in the entry region. This results in a flow that has
good cross sectional mixing, but minor longitudinal mixing, due to the

nearly uniform velocity profile Of a turbulent flow (Figure 14).

27

 

 

 

 

 

 

 

CS
H I5
E'
A
T
c
f»
A
c
/
7); /0
’meF
5
/
FIGURE 12.

26 30 32 36 40 44
TEMPEPA TURE °F

Soil Thermal Capacity as a Function Of Temperature

48

TurkfluqxxU Is

 

28

83 8:3 28.: 832:8 e3 .2 mesa:

     
 

GEES: $385: / dew

2MNQQK

 

 

 

     
  

20m.
2MNQQ...‘

205%: $355: 1\

\E «62me .92 \Nmmg.

 

.290 2MNQ~E
MQEQDM |\~

29

peso wee ea epeceea xeeeo_e> .e_ mmauea

 

m 2KOQQ
2.2925be

 

 

/1\

 

 

 

 

 

 

 

 

 

A\

sand «2

 

 

 

NOS a

form

is us

‘lalul

the

DISC

30

2.11 Numerical Stability
Finally, before the numerical model is developed the problem
was analyzed to determine if it would be numerically stable. The

formula
2 *
Ax /v At > 6

is used to insure three dimensional numerical stability using these

values

a = .00144 ftz/hr

Ax = 1 ft2

1/360 hours

At
the above criterion gives
Ax2/a*At = 250,000

predicting a numerically stable solution.

CHAPTER 3
NUMERICAL MODEL

The geometry and boundary conditions for the numerical model
are shown in Figure 15. A description of the procedures performed
in the numerical analysis, in the order in which they are performed
in the program follows, beginning with the values chosen for the

program variables.

3.1 Soil Properties and Initial Temperatures

Initially, reasonable values for the soil properties and initial
soil temperatures had to be determined (Table l). The deep earth
temperature in central lower Michigan is 50 degrees Fahrenheit (24).
The lowest average ambient temperature in central lower Michigan occurs
during the month Of January and is 30 degrees Fahrenheit (25). In
order to create a "worst case" condition for the input air temperature
to the duct, it was held at zero degrees Fahrenheit. The heat capacity
Of air is 0.24 BTU/1bm-°F and its thermal conductivity is 0.0142 BTU-l
HR-FT-°F (26). The overall coefficient of transmission from the soil
to the air is 1.40 BTU/HR-FT2-°F (27). The density of air at 20 degrees
Fahrenheit is 0.08275 1bm/FT3(28). Clay was chosen as the soil for the
simulation. The thermal conductivity of Healy clay with 10 and 20
percent moisture content is 0.458 and 0.833 BTU-FT/HR-FT2-°F, respectively,

and their respective heat capacities were .2 and .3 BTU/1bm-°F. The

31

32

Pete: Fcowcoszz

.m_ umame

mmo $2th

 

 

 

 

 

0: <m$Q<

 

MK

MSG 4m

x Q (Q23

Om

 

 

 

 

 

 

 

 

 

 

 

 

 

TL__

 

 

 

52.652 mmkaXT

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

m we; were L «$28 Mi
2 - 2
m. M.
2 2
m m
w .83: w.
a \ \ x
To. Q& q
a s
Q Q
2 z
a D
o o
m E
E2 m: o :2 8E
¢
N

Wk

 

$2qu

 

$52:

 

Em

 

 

 

 

 

 

 

 

Average
Current .
Deep 30‘;
Heat Ca:
Conduct‘
Overall
Density
NmnalE
Frozen E
Initial
Final MC
Frozen S
an13
3011 Der
Size Of
Size Of
Size of
Velocity
Number

Number 0
Number c
DGpth of
TemDerat
Number C
Latent L
Heat Ca;
Heat Ca;
—~..__~‘

 

33

Table 1

Values Used in the Finite-Difference Equations

 

Average Ambient Temperature
Current Ambient Temperature
Deep Soil Temperature

Heat Capacity Of Air
Conductivity of Air

Overall Heat Transfer Coefficient
Density of Air

Normal Soil Conductivity
Frozen Soil Conductivity
Initial Moisture Concentration
Final Moisture Concentration
Frozen Soil Heat Capacity
Normal Soil Heat Capacity
Soil Density

Size of X Increment

Size of Y Increment

Size of Z Increment

Velocity of Air in Pipe
Number of Y Increments
Number of Z Increments
Number of X Increments

Depth of Tube

Temperature Difference
Number of Time Steps

Latent Heat of Fusion

Heat Capacity of Water

Heat Capacity of Ice

TAVGAM
TCURAM
TSDILD

CPAIR =

KAIR
UAIR

KSOILN
KSOILF
IMOIST
FMOIST
CPSOLF
CPSOLN
DNSOIL
XINCRE
YINCRE
ZINCRE

N
Z
(Z
3
II II II II

TDIFFD

3
O
5
.2

O

O
4

.0142
1.40
DNAIR = .08275

II II II II II II II II II II II
m. C

5

I —J

No. a o

.5

O

0010'! U'INUJN—J

01

TIMEX = 432000

CPFUSE
CPNATR

CPICE =

1
1
.5

43.3
.0

BTU/LBM-°F 2
BTU-FT/HR-FT -°F
BTU/HR-FT2-°F
LBM/FT3
BTU-FT/HR-FT2-°F
BTU—FT/HR-FT2-°F
% by weight
% by weight
BTU/LBM-°F
BTU/LBM-°F
LBM/FT3

FT

FT

FT
FT/MIN

°F

BTU/LBM-°F
BTU/LBM-°F
BTU/LBM-°F

 

34

density of the clay is 95 lbm/FT3 (29). The heat capacity Of water
is 1.0 BTU/lbm-°F, the heat capacity of ice is .5 BTU/1bm-°F, and the
latent heat of fusion of water is 143.3 BTU/lbm-°F.

3.2 Soil Temperature Profile

After the thermal properties of the materials were determined, the
soil temperature array was initialized with a temperature distribution
in the z-direction. The initial temperature distribution is calculated

from the equation
Td-(Td-Ta)*(K/Z+1)2
where

Td is the deep earth temperature

Ta is the average ambient temperature

K is the distance from the deep earth isotherm

2 is the total depth being considered.
For a total depth of 11 feet, an average ambient temperature of 30
degrees Fahrenheit and a deep earth temperature of 50 degrees Fahrenheit

the temperature distribution is presented in Table 2. These temperatures

are loaded into the horizontal planes of the soil temperature array.

3.3 Determination of Time Step

The time step is determined by dividing the specified air velocity
by the element 1Ength in the pipe direction. Linking the time step to
the velocity and soil element size simplifies the convective term in
the heat transfer equation since the fluid is then convected one in-
crement per timestep and an insignificant amount of longitudinal mix-
ing can be assumed to occur during a time step due to the nearly uniform.

turbulent velocity profile.

35

Table 2

Initial Soil Temperature Distribution

 

Depth Temperature
(feet) (°F)

 

...a.—a.._a
NAOQmNO‘U'l-P-wN—‘O
U1
0‘
U1
(A)

 

36

3.4 Moisture Migration Function
A moisture attraction function, based solely on temperature,
will be used to account for increasing soil moisture with decreasing

temperature. The formula used is a second order polynomial Of the form
M = M. + (Mf-Mi)*(l-((TX-32)/(Td-32)2)
where

M is the moisture content Of the soil block being evaluated
Mi is the initial moisture content of the soil
Mf is the final moisture content Of the soil, and

T is the temperature of soil element.

3.5 Determination of Soil Heat Capacity and Conductivity

as a Function of Temperature

The next step in the program is the calculation of the thermal
properties of the soil as a function of temperature. There are two
important thermal properties considered, thermal conductivity and heat
capacity. Significant changes take place in these values in three
temperature ranges: above the freezing point Of water, below 30 degrees
Fahrenheit, and between 32 and 30 degrees Fahrenheit. In the program
the heat capacity of the soil, whose temperature is above freezing, is
calculated by summing the dry soil heat capacity and the heat capacity
of the moisture present in the soil. The heat capacity of frozen soil
is determined by summing the heat capacity of the dry soil and the
heat capacity of the ice present in the soil. The heat capacity of '
soil which is in the process of freezing is found by summing the heat

capacity of dry soil, the heat capacity of the water present in the soil,

37

and the latent heat of fusion Of the soil moisture (Figure 12).

The conductivity of the clay chosen for the simulation is .458
BTU/HR-FT-°F with 10 percent moisture content and .833 BTU/HR-FT-°F
with a 20 percent moisture content (30). Since the maximum moisture
content occurs at freezing and below, the 10 percent moisture content
conductivity was used down to 30 degrees Fahrenheit and the 20 percent
moisture content conductivity was used as the conductivity of frozen

soil below 30 degrees Fahrenheit.

3.6 Derivation of Heat Transfer Equations

There are three basic heat transfer equations which describe
the heat transfer in the model. They are: When the element is a
soil element surrounded by soil elements (see Figure 16) the follow-

ing equation is used:

- = = 9591. -
csms(Ti,j,k,2+l Ti,j,k,t) AE ks Ax (Ti-1,j,k,t Ti,j,k,£)

At + k éE—A-x-(T

Axay
s AX ) At + k (T.

i+l,j,k,t ‘ Ti,j,k,2 s AZ 1,j,k-1,2 ’

922x _ Ages.
Ti,j,k,£) At l ks AZ (Ti,j,k+i,t Ti,j,k,£) At T kx Ay
AXAZ
(Ti,j-i,k,t ' Ti,j,k,t) At + k5 Ay (Ti,j+i,k,t ' Ti,j,k,2) At

This is solved for T1 ,2+1 to give:

,J,k

38

mocmempm P_om an cmucaocczm mucmem_m Pwom op cmemcmch Hem:

.op mmzwmm

 

 

 

 

 

 

 

 

 

 

 

 

XAIIIII.

 

 

 

 

 

221.1 m
t
*

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

39

 

At kSAZAy
Ti,j,k,t+1 = csmx Ax (Ti-1,j,k,2 ' Ti,j,k,2) +
AZAy _ AAA! - -
ks Ax (Ti+l,j,k,2 Ti,j,k,t) ks AZ (Ti,j,k-i,t Ti,j,k,£) T

AXAy _ AXAz
ks Az (Ti,j,k+i,t Ti,i,k,t) T k (T

s Ay i,j-l,l,e ' Ti,j,k,2) T

AXAZ
ks7741 *1

thmz‘homa” LLLP

When the element borders on the pipe or the surface, such as an
element in contact with the top face of the tube (Figure 17), the follow-

ing equation is used:

- AE = k APAY-(T

CsmSTTi,0,l,2+l ’ Ti,0,l,2) ' s Ax i-l,0,l,2 ' Ti,0,l,2)

Asz _ -
At + k5 AX (T1+190919£ Tiaoslsl) At + UAXAY(TI9O’O’2 T1’0’1’l)

 

 

At T ks ATS—El (Ti,0,2,2 ' Ti,0,i,t) At T ks TAT? (TL-1,1,2 ' Ti,0,1,£)
At T ks TATE (Ti,l,1,2 ‘ Ti,O,l,t) At
Hhen solved for Ti 0,1,2+1 there results:
Ti,0,1,2+l T c:;s ks::Ay (Ti—1,0,1,t ' Ti,0,l,2) T ks TETT
(Ti+1,0,1,2 ' Ti,0,1,£) T UAXAy(11.0.0.2 ' Ti,0,1,2) T ks éfigx'
(Ti,0,2,t ' Ti,0,l,2) T ks T§$3'(Ti,-T,i,t ' Ti,O,1,£) T ks A232
‘(Hnnm‘TumLUJTHona

4O

HCMEmFm fiwzrm m mLmULom cow—ts HCQEm—m Pwom m :0 mucmpmm “mm: .N~ mmzwmm

 

 

 

 

 

 

 

 

 

 

 

 

 

\

 

.1 s m
w

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Nk—V-O:

41

When the element is the one that contains air inside the air tube,

the following equation is used (Figure 18)!

At -
camaua(Ti,O,O,t ‘ Ti-1,O,O,£) Ai'T Cama(Ti,O,O,e+l ’ Ti,0,0,2) ‘

= AZAy _ AZAy _
A5 ka Ax (Ti-1,0,0,2 Ti,O,O,£) At T ka Ax (Ti+1,0,0,t Ti,0,02)

At T UAXAY (11.0.1.2 ' Ti.0,0,£) At T UAXAV (Ti.0,-1.£ ' Ti,o,0,t)
At + UAXAz (Ti,1,0,2 - Ti’0,0,£)At -+ UAXAz (Ti,-l,0,t - Ti,0,0,2)

Solving Ti,0’0,£+1,_replacing At/Ax = l/ua, yields:

kaAZAy

 

= At _ AEAX
Ti,0,0,2+l Gama Ax (Ti-1.0.0.2 Ti.0.0,2) + k6 AX

(Ti+i,0,0,x ' 11.0.0.2) T ”AXAY (Ti,0,-l.2 ' Ti,0.0.2) +

UAXAY (Ti,O,1,t ' Ti,0,0,t) T UAXAZ (T ' T 1 T

i,-l,0,2 i,0,0,2

)] + T.

UAxAZ (Ti,],o,g ' Ti,o,0,e 1-1,0,0,2

3.7 Numerical Heat Transfer Equation
It should be noted that the finite difference equation in the
three modes is Of the same form and only the coefficients vary.
The basic finite difference equation used in the program to

calculate the heat transfer through an individual element is

42

 

 

 

 

 

 

 

 

 

 

 

 

Ne— v+¥

acmempm uwzpm a co mocepem Hem:

.w_ mmzomm

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

43

T(I,J,K,L+1) = KMAJOR*(K1*(T(I-1,J,K,L) - T(I,J,K,L)) +
K2*(T(I+1,J,K,L) - T(I,J,K,L)) + K3*(T(I,J,K-1,L) - T(I,J,K,L)) +
K4*(T(I,J,K+1,L) - T(I,J,K,L)) + K5*(T(I,J-1,K,L) - T(I,J,K,L)) +
K6*(T(I,J+1,K,L) - T(I,J,K,L))) + T(I,J,K,L).

The convection equation, used to describe the heat transfer to

the air in the pipe, is very simular
T(I,J,K,L+l) = KMAJOR*(K1*(T(I-1,J,K,L) - T(I,J,K,L)) +
K2*(T(I+1,J,K,L) - T(I,J,K,L)) + K3*(T(I,J,K-1,L) - T(I,J,K,L)) +
K4*(T(I,J,K+1,L) - T(I,J,K,L)) + K5*(T(I,J-1,K,L) - T(I,J,K,L)) +
K6*(T(I,J+1,K,L) - T(I,J,K,L))) + T(I-l,J,K,L)

with the only change in the last element reflecting the flow of air
in the tube.

This equation is used with different coefficients depending on
the position of the element. When the element is surrounded only by

soil elements, the coefficients used are

KMAJOR = AT/CS*MS

K1 = kS*AZ*AY/AX
K2 = k5*AZ*AY/AX
K3 = k5*AX*AY/AZ

44

K4 = kS*AX*AY/AZ
K5 = kS*AZ*AX/AY
K6 = kS*AZ*AX/AY

If an element is at the surface the 2+1 coefficient K4 is the air-

to-soil convection term
K4 = U*AX*AY

If the element is adjacent to the tube the y-l coefficient K5 must

be the air-to-soil convection term
K5 = U*AZ*AX

If the element is directly below the tube, the 2+1 coefficient K4 is

the air-to-soil convection term
K4 = U*AX*AY

If the element is directly above the tube, the 2-1 coefficient K3 is

the air-to-soil convection term
K3 = U*AX*AY

Finally, for the element containing the air inside the tube, the
heat capacity and mass used in the KMAJOR term are those of-air and

the conductivity used in the x+l and x-l terms is that of air:

KMAJOR = AT/Cf*Mf

Kl = ka*AZ*AY/AX

45

K2 = ka*AZ*AY/AX
K3 = U*AX*AY
K4 = U*AX*AY
K5 = U*AZ*AX
K6 = U*AZ*AX

3.8 Soil State Change Adjustment

After the heat transfer to an individual element is calculated,
the initial and final temperatures must be checked to determine if
the temperature passed through, into, or out of the 30 to 32 degree
Fahrenheit temperature range where the thermal properties of the
soil change (Figure 13). If this happens the temperature must be
adjusted in order to account for the changing thermal properties
since the temperature change calculated is based on the thermal
properties associated with the initial temperature. Because the
initial temperature could be in any one of the three ranges (below
30, above 32 and between 30 and 32) and the final temperature could
also be in any one of the three ranges there are nine possible com-
binations of initial and final temperatures.

In three of the cases, when the initial and final temperatures
are both in the same range, no adjustment is required.

In the other cases, the total heat transfer q, based on the
initial thermal properties, is calculated; then the amount of heat
required to bring the soil to the state change boundary is subtracted
from g, and finally the difference is divided by the thermal capacity

of the soil in the thermal range into which the temperature has passed.

46

This temperature change is added to or subtracted from the apprOpriate
boundary temperature. If the temperature has passed through the 30 to
32 degree range, the same adjustment is performed if the heat transfer

g calculated is less than the amount necessary to actually cause the
temperature to move through this range. Otherwise the final temperature

is calculated using the thermal properties of all three ranges (Table 3).

3.9 Preparation for Subsequent Time Step

After the heat transfer through the array has been calculated, a
number of adjustments are made before the next time step is begun. The
first is a test for output temperature stability. If the output
temperature is stable, the program exits the loop and notifies the user

that a stable heat transfer situation has been reached.

3.10 Adjustments for Insulation Condition and Plane of Symmetry

The insulated condition specified for the end planes must now be
established. The X=0 plane is filled with X=1 temperatures and
temperatures of the elements belonging to the plane just beyond the
end of the pipe (XNUM+1) are matched with those of the plane at the
end of the pipe (XNUM) (Figure 19). This simulates insulated planes
at each end of the pipe since no heat transfer will take place between
elements at the same temperatures. Similarly, the temperature Of the
elements on the negative (J-l) side of the pipe are filled with the
temperature of the element of the positive (J+1) side of the pipe
taking advantage of the plane of symmetry in the X-Z plane at J=0
(Figure 20). NO changes are made in the initial temperatures of the
elements more than 7.5 feet from the side of the tube or 7.5 feet below

the tube since the effects of the heat extraction are not expected to

 

k2w§<kmb§Q<
k2m2km5§Q< Q2 «TIE. m9 <2 2.92022an Nnmk

 

«TICK. mm2<Q .RZWEKWDEQY Q2 2+2: QMEQ.‘ NMVKVOM
. . zit

 

47

 

k2m§fimb§Q <
xm 402209 2 +2.2. 21$: Q .2 H2w2hwbfim< 02 on Wk

 

 

 

 

MWIATK NM Vk V Qm4 00m:
CK

a_aee Dessemane< aaeaeu aeaem

m m4m<b

48

eaeeeeeau ecu eaea_ameH ea eaeeeaem_eam .a_ “enema

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

" n ” \SQQK ‘1 . . _ _
2+E:QXNX.\\ \HX.\\
EBQXNX QNX
2.N.x.\+esext QM: 5L

1 Nstx‘EDQXVK. N HNfi-NVA‘QVF

49

 

 

 

 

 

 

 

 

 

 

I
J=-/ I J=I
l
I
1
12
,z”“‘L)LK:7-
T(X,-/,K,L+/) = +
T(X,/,K,L+I) Y
I
I
l
I
I
l
I
I
I
I
I
I
|
I
I
1
I
I
I
L__

PLANE OF S YMME TRY

FIGURE 20. Symmetrical Heat Transfer about the Plane
Y = 0

50

penetrate beyond 7.5 feet in the course of a heating

season.

3.11 Reinitialization of Soil Temperature Array
Finally, the t=0 array is filled with the newly calculated

t=l values and the next iteration is performed.

3.12 Programming Technique

Sections of the program were developed and tested independent
of the main program. The initial temperature distribution program
segment and thermal property state change adjustments were developed
as individual programs on a microcomputer to assure that they func-
tioned properly before being incorporated in the main program.

A large number of debugging statements were incorporated in
the code reporting on the progress of the calculations and the values
determined. Among these statement were reports on the initial
temperature distribution, which set of heat transfer coefficients
had been selected, and which of the nine possible initial-to-final
temperature adjustments had been selected.

After the program was complete, a version was written which
allowed single passes to be made through the program and initial
and final temperatures as well as position information to be inserted
to test the numerous nested if-blocks in the program in order to verify
the selection of the proper program segment. Using this version of
the program it was determined that the proper set of heat transfer
coefficients was selected for each position in the array and the
proper state change adjustment was selected for each possible combina-
tion of initial and final temperatures. A copy of this program is

included in Appendix B.

51

Finally, in order to speed the execution of the program the if-
blocks were arranged so that conditions were selected in order of
their frequency of occurrance. For instance, the most frequently
used set of heat transfer coefficients would be the soil-to-soil set
selected when the distance from the pipe is greater than or equal to
two increments in the y-direction. This is the first "if" encountered
after the surface elements are accounted for in the heat transfer
coefficient selection section of the program. In the interactive
version of the program error recovery was incorporated so that data
could be preserved even in the case of failures in "opens," "reads,"

"writes, and "closes" (Appendix C).

3.13 Computers Used

The heat transfer program was run on a CDC 750 and a Prime 7500.
A copy of the version used to produce the results can be found in
Appendix A. An interactive version for a microcomputer exists but
has not been tested. If a microcomputer could be obtained with a
floating point compressor and a Fortran compiler which used the co-
processor, the program could be run successfully on a microcomputer,

though 12 hour runs would probably be necessary.

52

CHAPTER 4
SUMMARY OF RESULTS

The simulation predicts that with zero degree Fahrenheit ambient
input air, the tube could be used to produce 240 cubic feet per minute
of 21.5 degree Fahrenheit air for one night, 120 cubic feet per minute
of 32 degree air for one night, or 60 cubic feet per minute of 40
degree Fahrenheit air overnight. The outlet temperature as a function
of time at 4 velocities is depicted in Figure 21. Figures 22, 23 and
24 show the temperature of the air versus position in the tube for various
times when the air velocities are one, two, and four feet per second,
respectively. The air temperature versus time at various x-positions in
the tube for air velocities of one, two and four feet per second are
shown in Figures 25, 26, and 27, respectively.

The computer model included the latent heat available during the
freezing of the soil; however, the sensible heat capacity of the soil
*was sufficient to provide the needed energy so that the first element
of soil never completely froze. Hence, an ice lense never formed

during any of the runs.

4.1 System Feasibility

The results of the simulation indicate that one duct with a
<:ross sectional area of one square foot could take in zero degree
Fahrenheit air and supply air at temperatures exceeding 40 degrees

Fahrenheit to the heat pump at the rate of 60 cubic feet per minute

53

 

 

 

 

mmwpvoopm> Lao; we we_h mamcm> ocsumcmaEmh uw—uzo .FN mmawmd
Amcnvp
~. __ o. m o u m m r n o
o.
Imam o.
:11 oN
mac 4
on
man N
ov
mac _

 

 

om

Aaav

54

vacuum can wood oco m_ xuvoopm> cw< on» cue:
weave msowce> com :o_upmoaix mzmcm> weaveLmQEmh .NN mmszm

 

Aummevx
om ow on ON 0. o
o.
\
.\
ON
on
.
e... a.
an...
a;

 

on H

55

accomm cog wood ozh mw xuwoopm> cw< on“ cog:
weave mzopcm> so; coeu_moaix mamcm> mczumcmasmh .mm mmzwfia

 

Apmmwvx
om ca on o~ o. 7 o

.v o.

o~

«I ~. on
m: a
x: v
m: _

ov

Aaav

 

on H.

56

ccouom cog Home caou we xu_oo_o> c_< asp cog:
woe?» maowce> com cora_moa-x mzmcm> mcaumchEmh .em umaoHd

 

Aummmvx
om ow on ON 0. o
0.
ON
:1 ~_
m: m
m: w
m: _
on
ov
CL
P

 

on

57

ccoomm com “com mco we xuwoopm> c_< esp cog: maze
one cw m:o_awmoa-x mzoece> no we?» mzmcm> acauecmasmh .mm mzswad

 

 

 

 

 

 

3.2:
N. : o. m o h m m t n N _ o
. 0.
Pump. 0.
ON
on
hum—u on I
9.
hum... om
Cay
A
om

58

uncomm son new; 03» m? xuwoo~m> sw< esp cog: maze
we“ cw mco_ppmoaix mzowce> an we?» mzmco> mezuecmasop .om mmzuHL

 

 

 

 

:2:
N. : o. m w h m w v. n N . o
o.
hum... 9
ON
kmm... on
om
hum... cm I
0?
CL
A
om

 

59

vacuum cog Home Lao; me xawoopm> cw< one can: mnap
msu cw meowuwmoaix maowce> on week mzmcm> mczaecoasmh .NN Manama

 

 

 

 

at:
N. : o. m o b m m f n N . o
hum: o.
o.
hum... on
ON
Hum... cm
on
O?
Gav
Om ._.

 

6D

for a 12 hour period. A 24,000 BTU heat pump would require approximately
1200 cubic feet per minute of 30 degree air to maintain its output (31).
Therefore, 20 ducts would be required to supply air to a heat pump of
this size. If the 50 foot long pipes were located 15 feet apart, assur-
ing no interference between pipes, 15,000 square feet of land would be
required. However, since sufficient heat is contained in only 2063
square feet of land to supply the heat pumps seasonal demand an area
substantially smaller than 15,000 square feet would probably be adequate.
If the full 15,000 square feet of land were required, an area 122 feet on

a side would be needed for pipe installation.

4.2 Economic Feasibility

The preceding analysis demonstrates that it is physically possible
to use a buried duct system to maintain the performance of the heat
pump in northern climates. However, the system must be economically
feasible. Hill the duct system save enough energy to pay for its
installation and upkeep? If a heat pump is supplied with 30 degree
Fahrenheit air, its coefficient of performance will be close to 2.8
(32). If a home requiring 83,000,000 BTU's of heat were heated with
electric resistance heat, 24,319 kilowatt hours of electricity would
be required. This electricity would cost $1824 at a cost of $.075
(33) per kilowatt hour. If a heat pump were used with an overall
coefficient of performance of 2.8, only 8685 kilowatt hours of electricity
‘would be needed at a cost of $651. This would yield an annual savings
of $1173.

Flexible irrigation pipe in the size required costs approximately

$2.00 per foot, so the pipe required would cost $2000. A heat pump

61

costs approximately $1000 more to install than central air conditioning
and electric resistance heating (34). A trench must be dug in which to
bury the pipe. Trenchers capable of digging 30 feet of trench per
hour rent for $30 per hour (35). Therefore, assuming free labor, the
cost of installing the duct would be $1000.

Finally, a small supplemental fan would be required to pull the
air through the duct system. The duct system has two 5 foot vertical
sections, a 50 foot horizontal section and an entry and an exit. Pulling
air through this system would create a pressure drop of approximately
0.02 inches of water (36). The fan required would consume about 50
watts of electricity and cost approximately $50 (37). During the course
of a heating season 333 kilowatt hours would be used, which at $.O75 a
kilowatt hour would cost $25.00. The system would cost $4050 more to
install than an electric resistance furnace and central air conditioning,
and would require a supplemental fan consuming $25 worth of electricity
per year. It would cost $1173 less to operate than electric resistance
heat each heating season, for a net annual savings of $1148. Assuming
an 18 percent interest rate the system would pay for itself in less

than 5 years.

CHAPTER 5
CONCLUSION

The program only allows the simulation of a duct system in the
single pass mode, and a very low constant ambient air temperature
was chosen for evaluation. If actual weather data were used and pro-
vision made in the program for Operation of the duct in a mode
appropriate for the current ambient temperature, system performance
would probably be significantly improved.

The buried duct system is an energy and cost effective method of
improving the performance of an air—to-air heat pump in a northern
climate. A five year payback is predicted when compared with an
electric resistance heating system. The economic analysis was per-
formed for a system using a duct air flow speed of 1 foot per second
and capable of supplying 100 percent of the heating requirement of a
residence at an ambient temperature of zero degrees Fahrenheit. Since
the duct system can supply 26 degree Fahrenheit air at 4 feet per
second and since the system need only operate down to ambient
temperatures of 10 degrees Fahrenheit in order to supply twice the
normal heat pump's contribution to the heating requirement, a much
smaller system could be installed which would yield an even higher

rate of return on investment.

62

CHAPTER 6
RECOMMENDATIONS

Since soil thermal properties for different soils vary greatly
and since different references give different values of thermal
properties for the same soil, experimental research should be carried
out to determine soil conductivity and thermal capacity as a function
of soil composition, soil temperature and soil moisture content. In
the program the soil conductivity was modeled as two discrete values
when in actuality the conductivity changes as the moisture content
of the soil changes and then increases by a factor of 1.4 when the
soil is completely frozen (38). Additionally, values as high as
1.6 BTU/FT-HR-°F were listed for clay with a 15 percent moisture content
(39).

Research should also be conducted to produce a curve fit which
will give the correct soil temperature distribution as a function of
the soil thermal properties, moisture content, ambient temperature and
time.

The process of moisture migration as a function of soil properties,
water table location, soil temperature and time should be studied and
a reliable estimator for this phenomenon should be determined.

The overall heat transfer coefficient from the atmosphere to the
surface should be determined for a variety of temperatures, wind speeds

solar insulation and surface conditions.

63

64

The overall heat transfer coefficient from the soil to the air
flowing in the pipe should be verified to determine the effects of
soil packing, pipe material and air velocity. Values as low as .4
and as high as 2 BTU/FTZ-HR-°F were noted in the literature. The
effect of pipe material is probably minimal since high density plastic
pipe has a conductivity of .226 BTU/FT-HR-°F and, in the .3 inch
thicknesses used, only affects the overall heat transfer coefficient
by about 3 percent when compared to higher conductivity materials such
as tile (40).

The program could be modified so that small elements could be
used in the vicinity of the pipe and large elements used away from
the pipe. This would speed the processing without significantly
affecting the accuracy of the simulation.

The number of cubic feet per minute of 30 degree Fahrenheit air
required for the heat pump's outside heat exchanger should be experi-
mentally verified.

A number of configurations could be simulated and tested such as
ducts with fins, rectangular ducts installed vertically and horizontally
with a number of different height-tO-width ratios (Figure 28).

Tests should also be made to find a formula for equating the
results for a rectangular duct to the heat transfer into a round pipe.
The program should be modified so that a weather record for a

typical winter complete with ambient temperatures, wind speeds,
humidity and solar insulation could be used in order to accurately
simulate an entire winter of operation. The program could be addi-

tionally modified to allow for Operation of the duct in single pass,

65

 

LS URFA CE

 

 

'“CONCRE TE BLOCK

 

 

 

 

a. Vertical Block Duct

 

 

 

 

 

 

b. Horizontal Block Duct

FIGURE 28. Possible Duct Configurations Using
Concrete Block

66

recirculating, bypass, and recovery modes depending on the current
weather conditions.

Modifications should be made to the program which allow it to
graphically display the temperature profile in the pipe as a function
of x-position and time step. Sufficient computer time should be
obtained so that a large number of runs could be made to optimize
the pipe length. Ducts of various sizes and shapes should be compared
in order to find the duct cross section which produces the most BTUs

per dollar.

APPENDICES

 

APPENDIX A
CYBER PROGRAM LISTING

67

1008

1208C
1308C
1408C
1508C
160'C
1708C
1803C
1908C
2003C
2103C
2203C
2308C
2408C
ZSOIC
2603C
2708C
2803C
290‘C
3003C
SIC-C
3203C
330=C
3403C
3508C
3603C
370-C
3803C
3908C
40080
410'C
420-C
4308C
4403C
4508C
4603C
470=C
460=C
4908C
500*C
510=C
5203C
5303C
5403C
5508C
5608C
570-C
5308C
5908C
6008C
610=C
6203C
630-C
6408C
6503C
6608C
670=C
6808C

68

PROGRAM TUBECYB
110-C§§*§§***§****§§*****a*§*§************§***********§***§§**§*§

PROGRAM TO DETERMINE THE THERMAL PERFORMANCE OF AN
UNDERGROUND TUBE

THIS PROGRAM PERFORMS A FINITE ELEMENT ANALYSIS

ON A BURIED TUBE THROUGH WHICH AIR IS FLOWING
CENTRAL DIFFERENCES ARE USED TO CALCULATE THE HEAT

GAIN OF AN ELEMENT AND FORWARD DIFFERENCES ARE USED

IN THE TIME STEP.

DATA FILES

TEMPFIL CONTAINS THE ARRAY TSOIL

VARIABLE LIST

TAVGAM AVERAGE AMBIENT TEMPERATURE

TCURAM CURRENT AMBIENT TEMPERATURE

TSOILD DEEP SOIL TEMPERATURE

CPAIR HEAT CAPACITY OF AIR

KAIR THERMAL CONDUCTIVITY OF THE AIR

UAIR SOIL-AIR HEAT TRANSFER COEFFICIENT

DNAIR DENSITY OF AIR

KSOIL THERMAL CONDUCTIVITY OF SOIL BASED ON ITS
CURRENT TEMPERATURE

KSOILN THERMAL CONDUCTRIVITY OF NORMAL SOIL

KSOILF THERMAL CONDUCTIVITY OF FROZEN SOIL

IMOIST INITIAL MOISTURE CONTENT OF SOIL

FMOIST FINAL MOISTURE CONTENT OF SOIL

CPSOIL HEAT CAPACITY OF THE SOIL BASED ON ITS
CURRENT TEMPERATURE

CPSOLF HEAT CAPACITY OF THE SOIL IN A FROZEN STATE

CPSOLN HEAT CAPACITY OF THE SOIL IN A THAHED STATE

CPSOLC HEAT CAPACITY OF THE SOIL BASED ON A COMBINATION
OF THANED HEAT CAPACITY AND HEAT OF FUSION

DNSOIL DENSITY OF THE SOIL

XINCRE THE INCREMENT SIZE IN THE X DIRECTION

YINCRE THE INCREMENT SIZE IN THE Y DIRECTION

ZINCRE THE INCREMENT SIZE IN THE Z DIRECTION

VELCTY THE VELOCITY OF THE AIR

XNUM NUMBER OF INCREMENTS IN THE X DIRECTION

YNUM NUMBER OF INCREMENTS IN THE Y DIRECTION

ZNUM NUMBER OF INCREMENTS IN THE Z DIRECTION .

TNUM THE DEPTH AT HHICH THE TUBE IS LOCATED VERTICALLY

TDIFFD TEMPERATURE DIFFERENCE REQUIRED FOR STABILITY

MOICON MOISTURE CONTENT OF SOIL AS A PERCENTAGE BY HEIGHT

TSOIL ARRAY CONTAINING TEMPERATURES OF SOIL BLOCK

TIMEX MAXIMUM NUMBER OF TIME STEPS TO BE PERFORMED

TEMPT TEMPERATURE OF PREPROCESSED INCREMENTS

SAVEAN ANSWER FOR DATA RETENTION

PREANS DIRECT CHOICE OF GENERATING INITIAL TEMPERATURES
USING THE PREPROCESSOR OR AN EXISTING TSOIL ARRAY

CPFUSE HEAT OF FUSION OF H20

CPWATR HEAT CAPACITY OF HATER

CPICE HEAT CAPACITY OF ICE

TIMSTP TIME STEP

MASSOL SOIL MASS DRY

MASSAR AIR MASS

MASS MASS

69

690-C KMAJOR COEFFICIENT OF HEAT TRANSFER EQUATION CONTAINING
700'C TIME STEP, HEAT CAPACITY AND SOIL MASS

710-C K1 COEFFICIENT OF THE X-1 TERM HEAT TRANSFER EQUATION
7208C K2 COEFFICIENT OF THE X+1 TERM HEAT TRANSFER EQUATION
7308C K3 COEFFICIENT OF THE 2-1 TERM HEAT TRANSFER EQUATION
7403C K4 COEFFICIENT OF THE 2+1 TERM HEAT TRANSFER EQUATION
750'C K5 COEFFICIENT OF THE Y-1 TERM HEAT TRANSFER EQUATION
760'C K6 COEFFICIENT OF THE Y+1 TERM HEAT TRANSFER EQUATION
7703C I LOOP COUNTER ARRAY INDEX

7BO‘C J LOOP COUNTER ARRAY INDEX

7903C K LOOP COUNTER ARRAY INDEX

8003C L LOOP COUNTER ARRAY INDEX

8108C M LOOP COUNTER ARRAY INDEX

3208C N LOOP COUNTER ARRAY INDEX

8303C T LOOP COUNTER ARRAY INDEX

B408C D LOOP COUNTER ARRAY INDEX

8503C B LOOP COUNTER ARRAY INDEX

8603C C LOOP COUNTER ARRAY INDEX

B70=C O LOOP COUNTER ARRAY INDEX

BBC-C TUBEHT TUBE LOCATION ABOVE 0 Z INCREMENT

8908C Q HEAT CAPACITY FOR STATE CHANGE CALCULATIONS

POO'C Q1 HEAT CAPACITY FOR STATE CHANGE CALCULATIONS

9103C Q2 HEAT CAPACITY FOR STATE CHANGE CALCULATIONS

9203C DELTAI TEMPERATURE CHANGE FOR STATE CHANGE CALCULATIONS
930-C DELTAZ TEMPERATURE CHANGE FOR STATE CHANGE CALCULATIONS
940-C DELTA3 TEMPERATURE CHANGE FOR STATE CHANGE CALCULATIONS
950‘C

9603C

970-C

9BO=C WRITTEN BY BRIAN LEARY

9908C AUGUST 10, 1984 VERSION

10008c§§§*§****§§*§§**§§**§*§§i§§§§§§§§§{‘{**§*§§*§§**§**§*§§*§§**§**
1010-C

1020=C SPECIFICATION STATEMENTS

1030-0

1040- IMPLICIT REAL (A-Z)

1050- integer i,J,k,1,m,n,xnum,ynum,znum,tnum,timex,
1060- +d,b,C,TIMLOP,D

1070- character savean

1080'II dimension tsoil(0:6,-1:8,0:12,0:1)
1090-C

1100-0 OPEN DATA FILE

1110-C

1120-20 open(6,2rr-4100,file-’TEMPFIL')
1130-C

1140-C INITIALIZE TEST VARIABLES

1130-C

1160-800 continue

1170- tavgam-SO

1180- tcuram-O

1190- tsoild=50 -
1200- cpair-.24

1210- kair=.0142

1220- uair-1.40

1230- dnair-.08275

1240- ksoiln-.4SB

1250- ksoilf-.833

1260- imoist-.1

1270- {moist-.2

1280- cpsolf-.3

1290a
13008
1310-
1320-
1330‘
1340'I
1350-
1360'
1370'I
1380-
1390-
14003
1410-
14203
1430-
1440a
1430-
1460-
14708
1480-
1490'
1500-1000
1310=C
1520=C
15308C
1540-
1550-1100
15603
1570-
1380a
1590-
1600-
1610-I
1620-
1630-1160
1640-1180
1650'I
1660=3
1670-
1680-
1690-1200
1700.
1710-C
1720=C
1730'C
1740-
17508C
1760=C
1770-C
1780-
1790=C
1800-C
1810-C
1820-
1830'C
1840'C
1850-C
1860'
1870-C
1880-C

7O

cpsoIn-.2
dnsoiI-95
xincre-lo
yincre-I
zincre-1
velcty-3600
ynum-7
znum-11
tnum-b
tdiffd-.0000000001
timex-36
TIHLOP-ioo
xnum-S
SAVEAN-'N'
preensai
cpfuse-143.3
cpuatr-1.0
cpice=.5
debug-0
debugz-O
I-O
continue

PREPROCESS ARRAY

if(preans.eq.2) then

read(*)tsoil

else

do 1200 k-0,znum+1
tempt-tsoild-((tsoild-tavgam)*((FLOAT(K)/FLOAT(ZNUH+1))**2))
do 1180 i=0,xnum+1

do 1160 j=-1,ynum+1

tsoil(i,j,k,1)=tempt
tsoil(i,j,k,1+1)-tempt

continue

continue

if(debug.ge.5) then

print *,'tempt= ',tempt,' at ',k

print *,'tsoi1- ’,tsoi1(1,0,k,1),' at ',k
endif

continue

endif

COHPUTE HEAT CAPACITY OF FREEZING SOIL
cpsolc=cpso1n+4moisticpuatr+4moist*cp¥use/2
COHPUTE POSITION OF TUBE IN SOIL
tubeht-znum-tnum+1

COHPUTE TIHE STEP

timstp-xincre/velcty

COMPUTE MASS OF AIR IN AN ELEMENT
massar-dnairfixincrefiyincreezincre

COMPUTE MASS OF SOIL IN AN ELEMENT

1390-0
1900-
1910-0
1920-0
1930-0
1940-
19:0-
1950-
1970-1800
1990-
1990-
2000-
2010-
2020-0
2030-0
2040-0
2050-2000
2060-
2070-
2080-
2090-
2100-
2110-0
2120-0
2130-0
2140-2010
2150-
2160-
2170-
2180-
2190-
2200-0
2210-0
2220-0
2230-0
2240-
2250-
2250-
2270-
2280-
2290-
2300-
2310-
2320-
2330-
2340-
2350-
2360-
2370-0
2380-C
2390-0
.2400-0
2410-0
2420-
2430-
2440-
2450-
2460-
2470-
2480-

 

71

massol-dnsoilfixincreiyincrefizincre
SET AIR TEMPERATURE AT TUBE ENTRANCE TO current ambient

do 1800 i-0,xnum+1
tsoil(i,0,tubeht,l)-(TCURAM)+((tavgam-tcuram)*
+((float(i)/float(xnum+1))**2))

continue

if(debug.ge.1) print *, 'preprocessing complete'
if(debug.ge.5.)print *, 'cpsoIc -',cpsolc,' tubeht =‘,tubeht,
+‘massar -',massar,' massol -‘,masso1

if<debug.ge.20)stop

LOOP CALCULATING HEAT TRANSFER THROUGHOUT THE ARRAY

continue

do 2400 t=1,timex
DO 2395 O=1,TIMLOP
do 2300 i=1,xnum
do 2200 j=0,ynum
do 2100 k=1,znum

CACULATE MOISTURE CONTENT AS A FUNCTION OF TEMPERATURE

i§(tsoil(i,j,k,1).gt.32)then
moicon-imoist+(fmoist-imoist)!
+(1-((tsoil(i,J,k,l)-32)/((tsoild-32)**2)))

else

moiconafmoist

endif ‘

SELECT SOIL HEAT CAPACITY AND THERMAL CONDUCTIVITY AS
A FUNCTION OF SOIL TEMPERATURE

if(tsoil(i,j,k,l).qe.32)then

cpsoil-cpsoln+moicon*cpwatr

ksoilaksoiln

else if(tsoil(i,J,k,1).le.30)then
cpsoil=cpso14+moiconicpice

ksoil-ksoilf

else

cpsoiI-cpsolc

ksoil-ksoiln

end if

if(debu92.ge.1) print 9, ‘soil thermal properties selected‘
if(debu92.ge.5)print *,'cpsoil -‘,cpsoil,'ksoil - ”,ksoil.
'moicon -',moicon

+

SELECT APPROPRIATE HEAT TRANSFER COEFFICIENTS
IF THE ELEMENT IS AT THE SURFACE USE AIR TO SOIL AT K+1

if(k.eq.ZNUM)then
kmaJor-timstp/(CPSOILfiMASSOL)
kl-ksoilizincrefiyincre/xincre
k2-ksoilizincreeyincre/xincre
kS-ksoilfixincreiyincre/zincre
k4-uairfixincre-yincre
k5=ksoi1*zincre*xincre/yincre

72

2490- kb-ksoiIfizincrefixincre/yincre

2500- if<debu92.ge.10) print *,‘choice - 1'

2510-:

2520-0 FOR ELEMENTS BEYOND THE FIRST TWO COLUMNS BUT BELOW
2530-0 THE SURFACE USE SOIL TO SOIL

2540-0

2550- else if(j.ge.2) then

2560- kmajor-timstp/(CPSOIL*MASSOL)

2570- k1-ksoi1*zincrefiyincre/xincre

2580- k2-ksoilezincreiyincre/xincre

2590- k3-ksoilfixincrefiyincre/zincre

2600- k4=ksoilfixincre§yincrel2incre

2610- ks-ksoilfizincre*xincre/yincre

2620- kb-ksoilfizincre*xincre/yincre

2630- i{(debu92.ge.10) print *,'choice - 2'

2640-0

2650-0 FOR ELEMENTS IN THE FIRST TWO COLUMNS BUT NOT CONTACTING
2660-0 THE TUBE USE SOIL TO SOIL COEF.

2670-0

2680- else if(k.ge.tubeht+2 .or. k.le.tubeht-2)then
2690- kmajor-timstp/(CPSOIL*MASSOL)

2700- kl-ksoilizincre*yincre/xincre

2710- k2-ksoil*zincre*yincre/xincre

2720- k3-ksoiIixincreeyincre/zincre

2730-

2740- k4-ksoil*xincre*yincre/zincre

2750- k5-ksoilfizincre*xincre/yincre

2760- kb-ksoilizincre*xincre/yincre

2770- if(debu92.ge.10) print *,'choice - 3'

2780-0

2790-0 ELEMENTS IN THE FIRST TWO COLUMNS ON A DIAGONAL FROM THE
2800-0 TUBE USE SOIL TO SOIL COEF.

2810-0

2820- else if(j.eq.1 .and. k.ne.tubeht) then

2830- kmajor-timstp/(CPSOILiMASSOL)

2840- kl-ksoilizincre*yincrelxincre

2850- k2=ksoilfizincrefiyincre/xincre

2860- k3-ksoi1*xincre*yincre/zincre

2870- k4-ksoil*xincre*yincre/zincre

2880- k5-ksoi1*zincreixincre/yincre

2890- kb-ksoilfizincre*xincrelyincre

2900- if(debu92.ge.10) print *,'choice - 4'

2910-0

2920-0 FOR ELEMENT ADJACENT TO THE TUBE USE AIR TO SOIL COEF
2930-0 FOR THE J-1 TERM

2940-0

2950- else if(j.eq.1 .and. k.eq.tubeht)then

2960- kmajor-timstp/(CPSOILQMASSOL)

2970- kl-ksoilezincre*yincre/xincre

2980- k2-ksoil*zincre-yincre/xincre

2990- k3-ksoilfixincre5yincre/2incre

3000- k4-ksoi1*xincre*yincre/zincre

3010- k5=uairizincreixincre

3020- kb-ksoilizincre§xincrelyincre

3030- if(debu92.qe.10) print *,'choice - 5'

3040-0

3050-0 FOR ELEMENT DIRECTLY BELOW THE TUBE USE AIR TO SOIL COEF.
3060-0 FOR THE K+1 TERM

3070-0

3080- else if(j.eq.0 .and. k.eq.tubeht-1)then

3090-
3100-
3110-
3120-
3130-
3140-
3150a
3160-
3170-0
3180-C
3190=C
3200-C
3210-
3220-
3230-
3240-
3250-
3260-
3270=
3280a
3290-
3300-C
3310=C
3320=C
3330-
3340=
3350-
3360-
3370-
3380-
3390-
3400=

3410-C
3420=C

3430'

3440-
3450-

3460-
3470-
3480-
3490s
3500-
3510-
3520-
3530-0
3540-C
3550-C

3560-2050

3570-
3580-
3590-
3600-
3610-
3620-
3630'C
3640-C
3650-0
3660-C
3670=C
3680=C

 

73

kmajor-timstp/(CPSOILiMASSOL)
kl-ksoilfizincre*yincrelxincre
k2-ksoilfizincre§yincre/xincre
k3-ksoilfixincreeyincre/zincre
k4=uair§xincrefiyincre
k5-ksoil*zincre*xincre/yincre
kbsksoilizincreexincre/yincre
if<debu92.qe.10) print *,'choice - 6'

FOR ELEMENTS DIRECTLY ABOVE THE TUBE USE AIR TO SOIL
FOR THE J-I TERM

else if(j.eq.0 .and. k.eq.tubeht+1)then
kmajor-timstp/(CPSOILfiMASSOL)
k1-ksoilfizincrefiyincre/xincre
k2=ksoilizincrefiyincre/xincre
k3=uair§xincrefiyincre
k4-ksoil*xincre*yincre/zincre
ks-ksoilizincreexincre/yincre
kb-ksoilezincre*xincre/yincre
if(debu92.ge.10) print *,'choice - 7'

COEFFICIENTS FOR THE ELEMENTS CONTAINING THE TUBE

else i§(j.eq.0 .and. k.eq.tubeht) then
kmajor-timstp/(CPAIR*MASSAR)
k1-kairfizincreeyincre/xincre
k2-kairfizincreiyincre/xincre
k3-uairexincre*yincre
k4=uairexincreeyincre
k5=uair§zincre§xincre
k6=uair§zincreixincre

CDNVECTIVE HEAT TRANSFER EQUATION

TSOIL(1,j,k,1+1)=kmajor*(k1*(TSOIL(i-1,j,k,I)-tsoil(i,j,k,1))+

+k2*(TSOIL(i+1,j,k,1)-tsoil(i,j,k,l))+
+k3*(TSOIL(i,j,k-1,l)-tsoil(i,J,k,l))+
+k4*(TSOIL(i,j,k+1,1)-tsoil(i,J,k,l))+
+k5*(TSOIL(i,J-l,k,1)-tsoil(i,j,k,l))+
+kb*(TSOIL(i,j+1,k,l)-tsoi1(i,j,k,l)))+TSOIL(I-1,J,K,L)

if<debu92.ge.10) print *,'choice - 8'

GO TO 2100

endif

if(debu92.ge.1) print *, 'heat transfer coefficients selected‘

HEAT TRANSFER EQUATION

TSOIL<i,j,k,1+1)-kmajor&(k1*(TSOIL(i-1,J,k,l)-tsoi1(i,j,k,l))+
+k2§<TSOIL(i+1,j,k,1)-tsoil(i,j,k,1))+
+k3§<TSOIL(i,j,k-1,1)-tsoil(i,j,k,1))+
+k4§<TSOIL(i,J,k+1,l)-tsoil(1,J,k,1))+
+k5fi<TSOILti,j-1,k,1)-tsoil(i,j,k,1))+
+k6i<TSOIL<i,J+1,k,l)-tsoil(i,j,k,1)))+TSOIL(I,J,K,L)

if(debu92.ge.1) print *,'heat trans€er equation completed'

ADJUST FOR TEMPERATURES CROSSING THE FREEZING ZONE OF
30 TO 32 DEGREES F

IF THE ELEMENT CONTAINS THE AIR PIPE NO ADJUSTMENT IS
NECESSARY

3690-0
3700-
3710-0
3720-0
3730-0
3740-
3750-
3760-:
3770-0
3780-0
3790-0
3800-
3810-
3820-
3830-0
3840-0
3850-0
3860-0
3870-
3880-
3890-
3900-
3910-
3920-
3930-
3940-
3950-
3960-
3970-
3980-
3990-
4000-0
4010-0
4020-0
4030-0
4040-
4050-
4060-
4070-
4080-
4090-
4100-0
4110-0
4120-0
4130-
4140-
4150-:
4160-0
4170-0
4180-0
4190-
4200-
4210-
4220-0
4230-0
4240-:
4250-0
4260-
4270-
4280-

74

if(j.eq.0 .and. k.eq.tubeht) go to 2100
IF THE INITIAL TEMPERATURE IS ABOVE 32

i+(tsoil(i,j,k,l).ge.32) then
if<debu92.ge.10) print *,‘l is greater than 32'

AND THE FINAL TEMPERATURE IS ABOVE 32 THEN NO
ADJUSTMENT IS NECESSARY

if(tsoil(i,j,k,l+1).ge.32) then
i§(debugz.ge.10) print *,'l+1 is greater than 32'
go to 2100

AND THE FINAL TEMPERATURE IS BELOW 30 THEN A
COMPOUND ADJUSMENT IS NECESSARY

else if<tsoil(i,j,k,l+1).le.30) then
i§(debu92.ge.10) print *,'l+1 is less than 30'
deltal-tsoil(i,j,k,l)-tsoil(i,j,k,l+1)
=cpsoln*delta1
q1=cpsoln*(tsoi1(i,j,k,l)-32)
q2=cpsolc§2

14(q-q1 .lt. q2)then
delta2=<q-q1)/:psolc
tsoil(i,J,k,l+1)=32-delta2

else

delta3-(q-q1-q2)/cpsoln
tsoil<i,j,k,l+1)-30-delta3

endif

IF THE FINAL TEMPERATURE~IS BETWEEN 30 AND 32 THEN THE
temperature is adjusted in THE FREEZING ZONE

else

if<debug2.ge.10) print *,'l+1 is between 30 and 32'
delta1-32-tsoil(i,J,k,l+1)
delta2-cpsoln/cpsolcfidelta1
tsoil(i,j,k,l+1)-32-delta2

endif

IF THE INITIAL TEMPERATURE IS BELOW 30 DEGREES

else i§(tsoil(i.j,k,l) .Ie. 30)then
if(debug2.ge.10) print *,'1 is less than 30'

AND THE FINAL TEMPERATURE IS BELOW 30 THEN NO ADJUSTMENT
IS NECESSARY

i{(tsoil(i,j,k,l+1).le.30)then
if(debug2.ge.10) print *,'l+1 is less than 30'
go to 2100

AND THE FINAL TEMPERATURE IS ABOVE 32 A COMPOUND adjustment
IS necessary

else if<tsoil(i.J,k,l+11.ge.32)then
if(debu92.ge.10) print *,'l+1 is greater than 32'
deltal-tsoil(i,J,k,l+1)-tsoil(i,j,k,l)

4290-
4300-
4310-
4320-
4330-
4340-
4350-
4360-
4370-
4380-
4390-0
4400-C
4410-C
4420-C
4430-C
4440-
4450-
4460-
4470-
4480-
4490-
4500=C
4510-C
4520-C
4530-
4540-
4550-0
4560=C
4570=C
4580-C
4590-
4600-.
4610-
4620-
4630-
4640-C
4650=C
4660-C
4670-C
4680-
4690-
4700-
4710-
4720-
4730-0
4740-C
4750-C
4760-C
4770-
4780-
4790-
4800-
4810-
4820-

4830-2100
4840-2200
4850-2300

4860-C
4870-C
4880=C

75

q-cpsolf*delta1
qi-cpsolf-(30-tsoil(i,j,k,l))
g2=cpsolc§2
i{(g-q1.1t.q2)then
deltaz-q-q1/cpsolc
tsoil(i,j,k,l+1)-30+delta2
else

delta3-(q-q1-q2)/cpsoln
tsoil(i,J,k,l+1)-32+delta3
endi+

AND THE FINAL TEMPERATURE IS BETWEEN 30 AND 32 THE FINAL
MUST BE REDUCED DUE TO THE HIGHER HEAT CAPACITY OF THE SOIL
IN THE FREEZING ZONE

else

if(debu92.ge.10) print *,'l+1 is between 30 and 32'
delta1-tsoil(i,j,k,l+1)-30
delta2-cpsolf/cpsolc*delta1
tsoil(i,J,k,l+1)-30+delta2

endif

IF THE INITIAL TEMPERATURE IS BETWEEN 30 AND 32

else
if(debu92.ge.10) print *,'l is between 30 and 32‘

AND THE FINAL TEMPERATURE IS ABOVE 32 DEGREES THE FINAL
TEMPERATURE IS INCREASED

if(tsoil(i,j,k,l+1).ge.32)then

if(debu92.ge.10) print *,'l+1 is greater than 32'
deltaI-tsoil(i,j,k,l+1)-32
delta2-cpsolc/cpsoln-delta1
tsoil(i,j,k,1+1)-32+delta2

AND THE FINAL TEMPERATURE IS BELOW 30 DEGREES THEN THE FINAL
TEMPERATURE MUST BE FURTHER REDUCED

else if(tsoil(i,j,k,l+1).le.30)then
if(debug2.ge.10) print *,'1+1 is less than 30'
delta1-30-tsoil(i,J,k,l+1)
deltaz-cpsolc/cpsolfideltal
tsoil(i,j,k,l+1)-30-delta2

AND THE FINAL TEMPERATURE IS BETWEEN 30 AND 32 DEGREES THEN
NO ADJUSTMENT IS NECESSARY

else

i+(debu92.ge.10) print *,'l+1 is between 30 and 32'
go to 2100 '

endif

endif

if(debu92.ge.1)print *,'state change accounted for'
continue

continue

continue

TEST FOR TEMPERATURE STABILITY

IF STABLE JUMP TO PROGRAM TERMINATION

4890-0
4900-
4910-
4920-
4930-
4940-0
4950-0
4960-0
4970-0
4980-0
4990-
5000-
5010-
5020-2310
5030-2320
5040-
5050-0
5060-0
5070-0
5080-0
5090-
5100-
5110-
5120-
5130-2330
5140-2340
5150-
5160-0
5170-0
5130-0
5190-
5200-0
5210-0
5220-0
5230-
5240-
250-
5250-
5270-2390
5290-
5290-2395
5300-
5310-
5320-
5330-
5340-
5350-2400
5360-0
5370-0
5380-0
5390-
5400-
5410-
5420-0
5430-0
5440-0
5450-2700
5460-2730
5470-
5480-0

76

if(t.gt.xnum+3)then
if(tsoil(xnum,0,tubeht,l+1)-tsoil(xnum,0,tubeht,l)
+ .lt. tdif+d) go to 2700

endif

TAKING ADVANTAGE OF THE PLANE OF SYMMETRY AROUND THE
TUBE IN THE X-Z PLANE FILL TSOIL(I,-1,K,L+l) ELEMENTS
WITH TSOIL(I,1,K,L+1) VALUES

do 2320 m-0,xnum+1

do 2310 n-0,znum+1
tsoil(m,-1,n,l+1)-tsoil(m,1,n,l+1)
continue

continue
if<debug.ge.1)print *, 'plane of symmetry filled'
MATCH END WALL TEMPERATURE TO ARRAY END TEMPERATURES
IN ORDER TO SATISFY THE INSULATION CONDITION

do 2340 J--1,ynum+1

do 2330 k-0,znum+1

tsoil(0,j,k,l+1)-tsoil(1,j,k,l+1)
tsoil(xnum+1,j,k,l+1)-tsoil(xnuM,J,K,L+1)

continue

continue

IF(DEBUG.GE.1.) PRINT *,'INSULATION CONDITION SATISFIED'

RESET TUBE ENTRANCE TEMPERATURE TO CURRENT AMBIENT
TSOIL(0,0,TUBEHT,L+1)=TCURAM
MOVE L+1 TEMPERATURES INTO THE L ARRAY FOR NEXT LOOP

do 2390 d-0,xnum+1

do 2390 b--1,ynum+1

do 2390 c-O,znum+1
tsoil(d,b,c,l)-tsoil(d,b,c,l+1)
continue
if(debug.ge.1)print *,
CONTINUE
hours-timstpfitimlop*t

'l+1 values inserted in l array'

print *, tsoil(xnum,0,tubeht,l+1), ' DEGREES AT '
+,HOURS,‘ HOURS.‘

WRITE(6,§) tsoil(xnum,0,tubeht,l+1), ‘ DEGREES AT '
+,HOURS,‘ HOURS.’

continue

IF MAXIMUM TIME STEP IS REACHED PRINT OUT MESSAGE

print i, ‘maximum time exceeded'

go to 2730

IF STABLE TEMPERATURE IS ACHIEVED PRINT MESSAGE
print 5, ‘output temperature stable at'

print 5, tsoil(xnum,0,tubeht,l+1), ‘degrees.’
print *, 'time- ',t

5490-C
5500-C
5510-
5520-
5530-
5540-3600
5550-
5560-
5570-C
5580-C
5590-C
5600-4100
5610-
5620-4300
5630-
5640-4500
5650-
5660-4700
5b70-
5680-

77

IF STORAGE WAS SPECIFIED WRITE OUT ARRAY VALUES

if(savean.eq.‘y‘)then
write(b,*,err-4500)tsoil

endif

:lose(b,err-4700)

print *, ‘normal termination'
STOP

ERROR MESSAGES

print *, 'open failed for data file'
STOP
print *,
STOP
print *,
STOP
print *,
STOP
END

‘read failed for data file'
'write failed for data file‘

'close *ailed for data file'

 

APPENDIX 8
TEST PROGRAM LISTING

[Ll

79

100- PROGRAM ONEPAS
110-0

120-0 SPECIFICATION STATEMENTS
130-0

140- IMPLICIT REAL (A-Z)

150- INTEGER I,J,K,L,M,N,XNUM,YNUM,ZNUM,TNUM,TIMEX,
160- +0.3,0

170- CHARACTER SAVEAN

180- DIMENSION TSOIL(0:b,-1:8,0:12,0:1)
190-0

200-0 OPEN DATA FILE

210-0

220-20 OPEN(6,ERR-4100,FILE-‘TEMPFIL‘)
230-0

240-0 INITIALIZE TEST VARIABLES
250-0

260-800 CONTINUE

270- TAVGAM-30

280- TCURAM-0

290- TSDILD-50

300- CPAIR-.24

310- KAIR-.0142

320- UAIR-1.40

330- DNAIR-.08275

340- KSOILN-.458

350- KSOILF-.833

360- IMOIST-.1

370- FMOIST-.2

380- CPSOLF-.3

390- CPSOLN-.2

400- DNSOIL-95

410- XINCRE-10

420- YINCRE-1

430- ZINCRE-1

440- VELCTY-3600

450- YNUM-7

460- ZNUM-11

470- TNUM-6

480- TDIFFD-.01

490- TIMEx-3

500- XNUM-5

510- SAVEAN-'N'

520- PREANs-1

530- CPFUSE-143.3

540- CPWATR-1.0

550- CPICE-.5

560- DEBUG-19

570- DEBU82-15

580- L-O

590-1000 CONTINUE

600-0

610-0 PREPROCESS ARRAY

620-0

630- IF(PREANS.EO.2) THEN
640-1100 READ<6,4)TSOIL

650- ELSE

660- DO 1200 K-0,ZNUM+1

670- TEMPT-TSOILD-((TSOILD-TAVGAM)*((FLOAT1K)/FLOAT(ZNUM+1))**2))

680- DO 1180 I-0,XNUM+1

690-
700-
710-1150
720-1180
730-
740-
750-
750-
770-1200
780-
790-0
800-0
810-0
820-
830-0
840-0
850-0
860-
870-0
880-0
890-0
900-
910-0
920-0
930-0
940-
950-0
960-0
970-0
980-
990-:
1000-0
1010-0
1020-
1030-
1040-
1050-
1060-
1070-0
1080-0
1090-0
1100-2000
1110-
1120-
1130-
1140-
1150-
1160-
1170-0
1180-0
1190-0
1200-2010
1210-
1220-
1230-
1240-
1250-
1260-0
1270-0
1280-0

80

DO 1160 J--1,YNUM+1
TSOIL(I,J,K,L)-TEMPT

CONTINUE

CONTINUE

IF(DEBUG.GE.5) THEN

PRINT -,'TEMPT- ‘,TEMPT,' AT ‘,K

PRINT -,'TSOIL- ‘,TSOIL(1,0,K,L),‘ AT ',K
ENDIF

CONTINUE

ENDIF

COMPUTE HEAT CAPACITY OF FREEZING SOIL
CPSOLc-CPSDLN+FMOISToCPNATR+FMOIST-CPFUSE/z
COMPUTE POSITION OF TUBE IN SOIL
TUBEHT-zNUM-TNUM+1
COMPUTE TIME STEP
TIMSTP-XINCRE/VELCTY
COMPUTE MASS OF AIR IN AN ELEMENT
MASSAR-DNAIRixINCRE§YINCRE§ZINCRE
COMPUTE MASS OF SOIL IN AN ELEMENT
MASSDL-DNSOIL-xINCRE*YIN0RE*21NCRE

SET AIR TEMPERATURE AT TUBE ENTRANCE TO CURRENT AMBIENT
TSOIL(0,0,TUBEHT,L)-TCURAM

IF(DEBUG.GE.1) PRINT 4, ‘PREPROCESSING COMPLETE‘
IF(DEBUG.GE.5.)PRINT 4, ‘CPSOLC -',0PSOLC,‘ TUBEHT -',TUBEHT,
+‘MASSAR -',MASSAR,' MASSOL -',MASSOL
IF(DEBUG.GE.20)STOP

LOOP CALCULATING HEAT TRANSFER THROUGHOUT THE ARRAY
CONTINUE

IF(DEBUG.GE.10) THEN

PRINT -,'INPUT I, a, K, L'

READ -, I,J,K,L

PRINT 4, ‘INPUT TSOILtL) AND TSOIL<L+1)'

READ 4, TSOIL(I,J,K,L),TSOIL(I,J,K,L+1)

ENDIF

CACULATE MOISTURE CONTENT AS A FUNCTION OF TEMPERATURE
IF(TSOIL(I,J,K,L).GT.32)THEN
MOICON-IMDIST+(FMOIST-IMOIST)4
+(1-<(TSOILII,J,K,L>-32)/<(TSDILD-3214-2111

ELSE

MOICON-FMOIST

ENDIF

SELECT SOIL HEAT CAPACITY AND THERMAL CONDUCTIVITY AS
A FUNCTION OF SOIL TEMPERATURE

81

1290-0

1300- IF(TSOILII,J,K,L).GE.32)THEN

1310- CPSOIL-CPSOLN4MOICON-CPNATR

1320- KSOIL-KSOILN

1330- ELSE IF(TSOIL(I,J,K,L).LE.30)THEN

1340- CPSOIL-CPSOLF+MOICON*CPICE

1350- KSOIL-KSOILF

1360- ELSE

1370- CPSDIL-CPSOLC

1380- KSOIL-KSOILN

1390- END IF

1400- IF(DEBU82.GE.1) PRINT 4, ‘SOIL THERMAL PROPERTIES SELECTED'
1410- IF(DEBUGZ.GE.5)PRINT 4,'0PSDIL -',CPSDIL,'KSOIL - ',KSOIL,
1420- +‘MOICON -',MOICON

1430-0

1440-0 SELECT APPROPRIATE HEAT TRANSFER COEFFICIENTS

1450-0

1460-0 IF THE ELEMENT IS AT THE SURFACE USE AIR TD SOIL AT K+1
1470-0

1480- IF(K.EO.ZNUM)THEN

1490- KMAJOR-TIMSTP/CPSOIL*MASSOL

1500- KI-KSOIL4zINCR84YINCRE/xINCRE

1510- K2-KSOIL*ZINCRE*YINCRE/XINCRE

1520- K3-KSOIL-xINCRE+YINCREIZINCRE

1530- K4-UAIR4xINCRE4YINCRE

1540- K5-KSOIL4zINCRE-XINCRE/YINCRE

1550- Ke-KSOIL-zINCRE4XINCRE/YINCRE

1560- IF(DEBU82.GE.10) PRINT -,'CHOICE - 1'

1570-0 FOR ELEMENTS BEYOND THE FIRST TWO COLUMNS BUT BELOW
1580-0 THE SURFACE USE SOIL TO SOIL

1590-0

1600- ELSE IF(J.GE.2) THEN

1610- KMAJOR-TIMSTP/CPSOIL*MASSOL
1620- Kl-KSOILillNCRE*YINCRE/XINCRE
1630- K2-KSOIL*ZINCREiYINCRE/XINCRE
1640- K3-KSOIL*XINCREiYINCRE/ZINCRE
1650- K4-KSOIL*XINCRE*YINCRE/ZINCRE
1660- K5-KSOIL*ZINCRE*XINCRE/YINCRE
1670- K6-KSOIL*ZINCRE*XINCRE/YINCRE
1680- IF(DEBUGZ.GE.10) PRINT *,'CHOICE - 2'
1690-0

1700-0 FOR ELEMENTS IN THE FIRST TWO COLUMNS BUT NOT CONTACTING
1710-0 THE TUBE USE SOIL TO SOIL COEF.

1720-0

1730- ELSE IF(K.GE.TUBEHT+2 .OR. K.LE.TUBEHT-2)THEN
1740- KMAJOR-TIMSTP/CPSOIL*MASSOL

1750- Kl-KSOIszINCRE§YINCREIXINCRE

1760- K2-KSOILQZINCRE*YINCRE/XINORE

1770- K3-KSOIL*XINCREQYINCRE/ZINCRE

1780-

1790- K4-KSOIL*XINCRE-YINCRE/ZINCRE

1800- K5-KSOIL*ZINCRE*XINCRE/YINCRE

1810- K6-KSOIL*ZINCRE*XINCRE/YINCRE

1820- IF(DEBU82.GE.10) PRINT *,'CHOICE - 3'
1830-0

1840-0 ELEMENTS IN THE FIRST TWO COLUMNS ON A DIAGONAL FROM THE
1850-0 TUBE USE SOIL TO SOIL COEF.

1860-0 .

1870- ELSE IF(J.EQ.1 .AND. K.NE.TUBEHT) THEN

1880- KMAJOR-TIMSTP/CPSOIL*MASSOL

1890- K1-KSOIL§ZINCRE*YINCRE/XINCRE

1900- K2-KSOILfiZINCREiYINCRE/XINCRE

1910- K3-KSOIL9XINCRE4YINCRE/ZINCRE

1920- K4-KSOILGXINCREfiYINCRE/ZINCRE

1930- K5-KSOILPZINCRE*XINCRE/YINCRE

1940- K6-KSOIL*ZINCRE*XINCRE/YINCRE

1950- IF(DEBU62.GE.10) PRINT *,'CHOICE - 4'
1960-0

1970-0 FOR ELEMENT ADJACENT TO THE TUBE USE AIR TO SOIL COEF
1980-0 FOR THE J-1 TERM

1990-0

2000- ELSE IF(J.EQ.1 .AND. K.EQ.TUBEHT)THEN
2010- KMAJOR-TIMSTP/CPSOILiMASSOL

2020- Kl-KSOIL*ZINCRE*YINCRE/XINCRE

2030- K2-KSOIL*ZINCRE*YINCRE/XINCRE

2040- K3=KSOIL*XINCRE*YINCRE/ZINCRE

2050- K4-KSOIL5XINCREiYINCRE/ZINCRE

2060- K5-UAIR*ZINCRE*XINCRE

2070- K6-KSOIL*ZINCRE*XINCRE/YINCRE

2080- IF(DEBUG2.GE.10) PRINT i,‘CHOICE - 5'
2090-0

2100-0 FOR ELEMENT DIRECTLY BELOW THE TUBE USE AIR TO SOIL COEF.
2110-0 FOR THE K*1 TERM

2120-0

2130- ELSE IF(J.EQ.0 .AND. K.EQ.TUBEHT-1)THEN
2140- KMAJOR-TIMSTP/CPSOILOMASSOL

2150- K1-KSOIL*ZINCRE*YINCRE/XINCRE

2160- K2-KSOILPZINCRE*YINCREIXINCRE

2170- K3-KSOIL*XINCRE*YINCRE/ZINCRE

2180- K4-UAIR*XINCRE*YINCRE

2190- K5-KSOILPZINCRE*XINCRE/YINCRE

2200- K6-KSOILiZINCREiXINCRE/YINCRE

2210- IF(DEBUG2.GE.10) PRINT *,'CHOICE - 6'
2220-0

2230-0 FOR ELEMENTS DIRECTLY ABOVE THE TUBE USE AIR TO SOIL
2240-0 FOR THE J-I TERM

2250-0

2260- ELSE IF(J.EQ.0 .AND. K.EQ.TUBEHT+1)THEN
2270- KMAJOR-TIMSTP/CPSOIL'MASSOL

2280- Kl-KSOIL*ZINCRE*YINCRE/XINCRE

2290- K2-KSOILGZINCRE*YINCRE/XINCRE

2300- K3=UAIR-XINCRE*YINCRE

2310- K4-KSOIL-XINCRE*YINORE/ZINCRE

2320- K5-KSOILGZINCRE*XINCRE/YINCRE

2330- K6-KSOILfiZINORE*XINCRE/YINCRE

2340- IF(DEBUGZ.GE.10) PRINT *,'0HOICE - 7'
2350-0

2360-0 COEFFICIENTS FOR THE ELEMENTS CONTAINING THE TUBE
2370-0

2380- ELSE IF(J.EQ.0 .AND. K.EQ.TUBEHT) THEN
2390- KMAJOR-TIMSTP/CPAIR*MASSAR

2400- Kl-KAIR*ZINCRE*YINCRE/XINCRE

2410- K2-KAIR-ZINCRE*YINCRE/XINCRE

2420- K3-UAIR5XINCREfYINCRE

2430- K4-UAIR*XINCRE#YINCRE

2440- K5-UAIR*ZINCRE*XINCRE

2450- K6-UAIR*ZINCRE*XINCRE

2460-0 CONVECTIVE TRANSFER EQUATION

2470-0

2480-2050 TSOIL(I,J,K,L+1)-KMAJOR*(K1*(TSOIL(I-l,J,K,L)+TSOIL(I,J,K,L))-

82

2600-0

2610-2050

2620-
2630-
2640-
2650-
2660-
2670-
2680-
2690-
2700-
2710-
2720-
2730-
2740-0
2750-0
2760-0
2770-0
2780-0
2790-0
2800-0
2810-
2820-0
2830-0
2840-0
2850-
2860-
2870-0
2880-0
2890-0
2900-0
2910-
2920-
2930-
2940-0
2950-0
2960-0
2970-0
2980-
2990-
3000-
3010-
3020-
3030-
3040-
3050-
3060-
3070-
3080-

83

+K2-<TSOIL(I+1,J,K,L>+TSOIL(I,J,K,L)1+
+K3§<TSOIL(I,J,K-1,L)+TSOIL(I,J,K,L)1+
+K44<TSDIL(I,J,K+1,L)+TSOIL(I,J,K,L)1+
+K54<TSOIL(I,J-1,K,L)+TSOIL<I,J,K,L)1+
+K54<TSDIL<I,J+1,K,L>+TSOIL(I,J,K,L1I)+TSOIL(I—1,J,K,L)

IF(DEBU62.GE.10) PRINT -,'CHOI0E - 8'

GO TO 2100

ENDIF

IF(DEBUGZ.GE.1) PRINT 4, ‘HEAT TRANSFER COEFFICIENTS SELECTED

HEAT TRANSFER EQUATION

TSOIL<I,J,K,L+1)-KHAJOR-(KI-(TSDIL1I-1,J,K,L1+TSOIL<I,J,K,L>>-
+K2-(TSDILII+1,J,K,L)+TSOIL(I,J,K,L)1+
+K3-(TSDILII,J,K-1,L)+TSDIL(I,J,K,L)1+
+K44<TSOIL<I,J,K+1,L)+TSOIL<I,J,K,L)1+
+K54<TSOIL<I,J-1,K,L)+TSOIL(I,J,K.L11+
+K6*(TSOIL(I,J+1,K,L)+TSOIL(I,J,K,L)))+TSOIL(I,J,K,L)

IF(DEBU82.GE.1) PRINT *,‘HEAT TRANSFER EQUATION COMPLETED‘

IFmamEGEIm flfiN

PRINT *,'INPUT 1, J, K, L‘

READ 4, I,J,K,L

PRINT 4, 'INPUT TSOIL(L) AND TSOIL(L+1)‘

READ 4, TSOIL<I ,J,K,L) ,TSOIL(I.J,K,L+1)

ENDIF

ADJUST FOR TEMPERATURES CROSSING THE FREEZING ZONE OF
30 TO 32 DEGREES F

IF THE ELEMENT CONTAINS THE AIR PIPE NO ADJUSTMENT IS
NECESSARY

IF(J.EQ.0 .AND. K.EQ.TUBEHT) GO TO 2100
IF THE INITIAL TEMPERATURE IS ABOVE 32

IF(TSOIL(I,J,K,L).GE.32) THEN
IF(DEBUG2.GE.10) PRINT *,'L IS GREATER THAN 32'

AND THE FINAL TEMPERATURE IS ABOVE 32 THEN NO
ADJUSTMENT IS NECESSARY

IF(TSOIL(I,J,K,L+1).GE.32) THEN
IF(DEBUG2.GE.10) PRINT §,'L+1 IS GREATER THAN 32'
GO TO 2100

AND THE FINAL TEMPERATURE IS BELOW 30 THEN A
COMPOUND ADJUSMENT IS NECESSARY

ELSE IF(TSOIL(I,J,K,L+1).LE.30) THEN
IF(DEBUG2.GE.10) PRINT -,'L+1 IS LESS THAN 30'
DELTAl-TSOIL(I,J,K,L)-TSOIL(I,J,K.L+1)
D-CPSDLN4DELTA1
Ol-CPSOLN*(TSOIL(I,J,K,L)-32)
Dz-CPSOL042

IF(Q-Ql .LT. 82>THEN
DELTA2-(Q-O1)/0PSOLC
TSOIL(I,J,K,L+1)-32-DELTA2

ELSE

DELTA3-(O-Ol-02)/CPSOLN

3090-
3100-
3110-0
3120-0
3130-0
3140-0
3150-
3160-
3170-
3180-
3190-
3200-
3210-0
3220-0
3230-0
3240-
3250-
3260-0
3270-0
3280-0
3290-
3300-
3310-
3320-0
3330-0
3340-0
3350-
3360-
3370-
3380-
3390-
3400-
3410-
3420-
3430-
3440-
3450-
3460-
3470-
3480-0
3490-0
3500-0
3510-0
3520-
3530-
3540-
3550-
3560-
3570-
3580-0
3590-0
3600-0
3610-
3620-
3630-0
3640-0
3650-0
3660-
3670-
3680-

84

TSOIL<I,J,K,L+1)=3o—DELTA3
ENDIF

IF THE FINAL TEMPERATURE IS BETWEEN 30 AND 32 THEN THE
IN THE FREEZING ZONE

ELSE

IF(DEBUG2.GE.10) PRINT *,'L+1 IS BETWEEN 30 AND 32'
DELTA1-32-TSOIL(I,J,K,L+l)
DELTA2-0PSOLN/CPSOLC*DELTA1
TSOIL(I,J,K,L+1)-32-DELTA2

ENDIF

IF THE INITIAL TEMPERATURE IS BELOW 30 DEGREES

ELSE IF(TSOIL(I,J,K,L) .LE. 30)THEN

IF(DEBUG2.GE.10) PRINT *,'L IS LESS THAN 30'

AND THE FINAL TEMPERATURE IS BELOW 30 THEN NO ADJUSTMENT
IS NECESSARY ~

IF(TSOIL(I,J,K,L+1).LE.30)THEN
IF(DEBUG2.GE.10) PRINT *,'L+1 IS LESS THAN 30'
GO TO 2100

AND THE FINAL TEMPERATURE IS ABOVE 32 A COMPOUND IS NECESSARY

ELSE IF(TSOIL(I,J,K,L+1).GE.32)THEN

IF(DEBUG2.GE.10) PRINT *,'L+l IS GREATER THAN 32'
DELTAl-TSOIL(I,J,K,L+1)-TSOIL(I,J,K,L)

Q-CPSOLFfiDELTAl _

Ql-CPSOLF*(30-TSOIL(I,J,K,L))

Q2-CPSOLC*2

IF(Q-Ql.LT.Q2)THEN

DELTA2-Q-Q1/CPSOLC

TSOIL(I,J,K,L+1)-30+DELTA2

ELSE

DELTA3-(Q-Ql-Q2)/CPSOLN

TSOILtI,J,K,L+1)-32+DELTA3

ENDIF

AND THE FINAL TEMPERATURE IS BETWEEN 30 AND 32 THE FINAL
MUST BE REDUCED DUE TO THE HIGHER HEAT CAPACITY OF THE SOIL
IN THE FREEZING ZONE

'ELSE

IF(DEBU82.GE.10) PRINT -,'L+1 IS BETWEEN 30 AND 32'
DELTAl-TSOIL(I,J,K,L+1)-30
DELTAz-CPSDLF/CPSOL0-DELTA1
TSOIL(I,J,K,L+1)-30+DELTA2

ENDIF

IF THE INITIAL TEMPERATURE IS BETWEEN 30 AND 32

ELSE

IF(DEBUGZ.GE.10) PRINT *,'L IS BETWEEN 30 AND 32'

AND THE FINAL TEMPERATURE IS ABOVE 32 DEGREES THE FINAL
TEMPERATURE IS INCREASED

IF(TSOIL(I,J,K,L+1).GE.32)THEN
IF(DEBUG2.GE.10) PRINT *,'L+1 IS GREATER THAN 32'
DELTAl-TSOIL(I,J,K,L+1)-32

3690-
3700-
3710-0
3720-0
3730-0
3740-0
3750-
3760-
3770-
3780-
3790-
3800-0
3810-0
3820-0
3830-0
3840-
3850-
3860-
3870-
3880-

3890-2100

3900-
3910-
3920-
3930-
3940-
3950-
3960-
3970-
3980-0
3990-0
4000-0
4010-0
4020-
4030-
4040-0
4050-0
4060-0
4070-0
4080-0
4090-
4100-
4110-

4120-2310
4130-2320

4140-
4150-0
4160-0
4170-0
4180-0
4190-
4200-
4210-
4220-

4230-2330
4240-2340

4250-

4260-0
4270-0
4280-0

85

DELTA2-0PSOLCICPSOLN*DELTA1
TSOIL(I,J,K,L+l)-32+DELTA2

AND THE FINAL TEMPERATURE IS BELOW 30 DEGREES THEN THE FINAL
TEMPERATURE MUST BE FURTHER REDUCED

ELSE IF(TSOIL(I,J,K,L+1).LE.30)THEN
IF(DEBUG2.GE.10) PRINT 4,‘L+l IS LESS THAN 30'
DELTA1-30-TSOIL(I,J,K,L+l)
DELTA2-0PSOLC/CPSOLF-DELTA1
TSOIL(I,J,K,L+1)-30-DELTA2

AND THE FINAL TEMPERATURE IS BETWEEN 30 AND 32 DEGREES THEN
NO ADJUSTMENT IS NECESSARY

ELSE

IF(DEBUGZ.GE.10) PRINT -,'L+1 IS BETWEEN 30 AND 32'
GO TO 2100

ENDIF

ENDIF

CONTINUE

IF(DEBU82.GE.1)PRINT 4,'STATE CHANGE ACCOUNTED FOR‘
IF(DEBUG2.GE.15)THEN

PRINT -, ‘ENTER 1 TO TEST ANOTHER CONDITION‘

PRINT -, 'ENTER 2 TO TERMINATE'

READ 4, LOOP

IF(LOOP.EQ.1) GO TO 2000

ENDIF

PRINT 4, TSOIL(XNUM,0,TUBEHT,L+1), ‘DEGREES.’

TEST FOR TEMPERATURE STABILITY
IF STABLE JUMP TO PROGRAM TERMINATION

IF(TSOIL(XNUM,0,TUBEHT,L+1)-TSOIL(XNUM,O,TUBEHT,L)
+ .LT. TDIFFD) GO TO 2700

TAKING ADVANTAGE OF THE PLANE OF SYMMETRY AROUND THE
TUBE IN THE X-Z PLANE FILL TSOIL(I,-1,K,L+l) ELEMENTS
WITH TSOIL(I,1,K,L+1) VALUES

DO 2320 M-1,XNUM

DO 2310 N-1,ZNUM
TSOIL(M,-1,N,L+1)-TSOIL(M,1,N,L+1)

CONTINUE

CONTINUE

IF(DEBUG.GE.1)PRINT *, 'PLANE OF SYMMETRY FILLED'

MATCH END WALL TEMPERATURE TO ARRAY END TEMPERATURES
IN ORDER TO SATISFY THE INSULATION CONDITION

DO 2340 J-1,YNUM
DO 2330 K-1,ZNUM

TSOIL(O,J,K,L+1)-TSOIL(1,J,K,L+1)
TSOIL(XNUM+1,J,K,L+1)-TSOIL(XNUM,J,K,L+1)

CONTINUE

CONTINUE

IF(DEBUG.GE.1) PRINT -,'INSULATIION CONDITION SATISFIED‘

RESET TUBE ENTRANCE TEMPERATURE TO CURRENT AMBIENT

86

4290- TSOIL(0,0,TUBEHT,L+1)=TOURAM
4300- IF(DEBUG.GE.1)PRINT 4,'AIR ADVANCED IN TUBE‘
4310-0

4320-0 MOVE L+1 TEMPERATURES INTO THE L ARRAY FOR NEXT LOOP
4330-0

4340- DO 2390 D-o,XNUM+1

4350- DO 2390 B--1,YNUM

4360- DO 2390 0-1,ZNUM

4370- TSOIL(D,B,C,L)-TSOIL(D,B,0,L+1)

4380-2390 CONTINUE

4390- IF(DEBUG.GE.1)PRINT 4, 'L+1 VALUES INSERTED IN L ARRAY‘
4400-2400 CONTINUE

4410-0

4420-: IF MAXIMUM TIME STEP IS REACHED PRINT OUT MESSAGE
4430-0

4440- PRINT 4, ‘MAXIMUM TIME EXCEEDED’

4450-

4460- GO TO 2730

4470-0

4480-: IF STABLE TEMPERATURE IS ACHIEVED PRINT MESSAGE
4490-0

4500-2700 PRINT 4, ‘OUTPUT TEMPERATURE STABLE AT'
4510-2730 PRINT 4, TSOIL(XNUM,0,TUBEHT,L+1), ‘DEGREES.’
4520- PRINT 4, 'TIMEx- ‘,TIMEX

4530-0

4540-0 IF STORAGE WAS SPECIFIED WRITE DUT ARRAY VALUES
4550-0

4560- IF(SAVEAN.EQ.'Y‘)THEN

4570- NRITE<6,4,ERR-45001TSDIL

4580- ENDIF

4590-3600 CLOSE(6,ERR-4700)

4600- PRINT 4, ‘NORMAL TERMINATION'

4610- STOP

4620-0

4630-0 ERROR MESSAGES

4640-0

4650-4100 PRINT 4, 'OPEN FAILED FOR DATA FILE'
4660- STOP

4670-4300 PRINT 4, 'READ FAILED FOR DATA FILE‘
4680- STOP

4690-4500 PRINT 4, ‘WRITE FAILED FOR DATA FILE'
4700- STOP

4710-4700 PRINT 4, ‘CLOSE FAILED FOR DATA FILE'
4720- STOP

4730- END

APPENDIX C
MICROCOMPUTER PROGRAM LISTING

87

88

*¥#*#######**#####¥##$##$*#$########$§$###*###*##N#$####
Program to determine the thermal Performance of an
underground tube
This Program Performs a finite element analysis
on a buried tube through which air is flowing
Central differences are used to calculate the heat

fluirlflclrsnlfirnflrlnrn

fir]nifirlnflflfintfifnflrlnii

fl '

Flnlfirlnlifnu

HHHH

1f1filfif|flf1flflflflfllfi
(VCJH kiH

‘Fafi

Gain of an element and forward differences

in the

are used
time steP.

Data Files

TemPfil
Paramfl

Variable List I Integer

HNSNEP
TRVBRM
TCURAM
TSDILD
CPAIR
KAIR
UAIR
DNAIR
KSOIL

KSOILN
KSOILF
IMOIST
FMOIST
CPSOIL

CPSOLF
CPSOLN
CPSULC

DNSOIL
XINCRE
YINCRE
ZINCRE
VELCTY
KHUM
TNUM
ZNUM
TNUM

TDIFFD
MOICON
TSOIL
TIMER

TCHRHG

SELECT
SRVERH
PRRRH

PRERHS

CPFUSE

SPURTR
CPICE

Contains the array tsoil
Contains soil thermal characteristic
dimension Parameters

and

C Char

V or H check on input Parameters
Rveraee ambient temPerature
Current ambient temPerature

DeeP soil temPerature

Heat caPacits of air

Thermal conductivity of the air
Soil~air heat transfer coefficient
Densits of air

Thermal conductiuite of soil based
current temPeraturefiD) '
Thermal conductivits of normal soil
Thermal conductivity of frozen soil
Initial moisture content of soil
Final moisture content of soil
Heat caPacits of the soil based
current temPeratureCCb
Heat caPacits of the soil
Heat caPacits of the soil in a thawed state

Heat capacity of the soil based on combinationCC)
of thawed heat caPacite and heat of fusion
Densits of the soil

The increment size in the x direction

The increment size desired in the Y direction
The increment size desired in the z direction
The velocity of the air
Number of increments in the
The number of increments in
The number of increments in
The dePth at which the tube
is located verticalla
TemPerature difference refluired for stability
Moisture content of soil as a Percentage by weightIC)
Rrras containing temPeratures of soil block

Maximum number of time stePs to be Performed
TemPerature change between PreProcessed incrementst?)
Menu choice for data inPut

Answer for data retention

Rnswer for Parameter retention

Direct choice of Generating initial temPeratures
using the PreProcessor or an existing tsoil arras
Heat of fusion of H20

Heat caPacitS of water

Heat caPacita of ice

on its

on its

in a frozen state

X direction
the Y direction
the 2 direction

Flfllifintln

tanfin

flfiifirln ni1flnififlfllfifint1€1flnlirlnu1rann1rb

lrnn

f'

HHHHHHHHHHH

TIMSTP
MRSSOL
MHSSHR
MESS

KMHJDR

C
TUBEHT
E

Q1

Q2
DELTHI
DELTHZ
DELTHB
H
TEMPT

DEBUG

DEBUGZ

Written be Brian

Time
Soil

stePCC)
mass

Rir massCC)
MassiC)E?]
Coefficient of

timestep,
Coefficient
Coefficient
Coefficient
Coefficient
Coefficient
Coefficient

LooP
LOOP
LooP
Loop
LooP
Loop
Loop
Loop
LooP
Loop
Tube
Heat
Heat
Heat

TemPerature change for state change
TemPerature change for state change
state change

heat
of
of
of
of
I) ‘F
of

IZif‘S ': I:

89

J

heat transfer enuation containin?

capacitai

the
the
the
the
the
the

K-1
2+1
2-1
2+1
v~1
V—i

term
term
term
term
term
term

counter;
counter;
counter,
counter;
counter;
counter;
counter,
counter,
counter;
counter,
location
capacita
caPacita
capacity

arras
array
arras
arraa
arras
array
array
arras
arras
arras
above

for state

for s
for 5

Temperature change
Hnswer for redoing i/o
Variable into which initial temperature
distribution is written

Debugging
different
Debu99in9
different

Modification History

HUBUST 7;

variable
levels of debugging output

Leary

1984 version

index
index
index
index
index
index
index
index
index
index
8 Z

tate
tate

for

and soil

heat
heat
heat
heat
heat
heat

increment

change calculations
change calculations
change calculations
calculations
calculations
calculations

mass

transfer
transfer
transfer
transfer
transfer
transfer

eQuatiOn
enuation
eQuation
eQuation
eQuation
equation

variable different values generate
levels of debugging outPut

different values 9enerate

caasaawawaeewaaaaeeaweaa##ewaeweaeeeaaeaeaeaeeaaeaaaaaaaa
SPecification Statements

*0
G

ntinufiflnlo
(3

IMPLICIT REEL

cn-zn

INTEGER I;J;K»L;MaN;KNUN;VNUMJZNUH;TNUMJTIMER,
SELECTxTiniBJC

CHHRHCTER

SRVEHHJHHSNERJHJPHRHH

DIMENSION TSDILCB=51i-132138331:B l)

DPen InPut data files

DPENCS;ERR=4BBBaFILEaPRRHMFLP
DPENCS;ERR=4IBB»FILE=TEMPFIL?

InPut

Determine inPut

flsk

user to choose

section

mode

whether

brifiufi

Q

001

unlfi

“[10

(Jfli1011fl

G

nrnn nrin

'fl

h

90

distribution data is to be read From a File or
numerically generated

HRITEifib ’Do eou wish to’
NEITEfifl) ’ ’

MRITE(#) ’1 Produce an initial temperature distribution’
HRITECW) ’numericalle?’

HRITE<¥> ’2 Read the data From a data File?’
NRITE(#), ’9.9 Terminate Program?’

WRITEC#). ’ ’

MRITEfifi); ’InPut the number c0rresPondin9 to aour’
MRITEcfi), ’selection’

REHDC$> PRERNS

IF option 1 or 2 is chosen then go to 56

IF(PREHH3.EQ.1) THEN

60 TD 58
ELSE IF(PREHHS.EQ.2) THEN
GO TD 58

I? termination of program is chosen then end

ELSE IFfiPRERHS=999> THEN
END

allow user to correct For keyboard input errors

ELSE

HRITECfi); ’PLEHSE SELECT R NUMBER FROM THE MENU’
GO TU 48

ENDIF

fisk user to choose whether the parameters are
to be read From the test data set, From a File or
input From the keyboard

NRITEfifi) ’Do sou wish to’

HRITECfi) ’ ’

MRITE(#) ’1 Read From the set of test data?’
NRITE(#) ’2 Read the data From a data File?”
NR1TEC$D ’3 Input the data via the keeboard?’
HRITEfifi); ’999 Terminate”

NRITEflfi), ’ ’

WRITEfifl), ’InPut the number corresPondin9 to aour’
HRITEfifi); ’selection’

RccePt menu choice €rom keyboard
REHDCfi) SELECT
I? test File is chosen move to line number 886

IF(SELECTSI) THEN
00 T0 886

I? data tile is chosen read Parameters

ELSE IF(SELECT=2) THEN

GI

I'lfln ++++“~i

fll'l n

n

HOOD

91

REHDfi5,ERRa4286;EunslooobTRVBHm,TCURHM,TSDILD,CPHIR,
KHIRTURIR)DHHIRJKSDILHJKSDILF,IMOISTJFMDIST,CPSDLFJ
CPSOLN,DNSOILax1HCRE,VIHCRE,ZINCREJWELCTVJVNUM,ZHUM,
THUMJTDIFFanIMEH,DSOILH;DEDILF,KHUHJSRWEHH,PRERHS,
CPFUSE,CPHRTR,DPICE

50 TD 1838

IF keyboard input is specified then move to line 188

ELSE IF(SELECT=3) THEN
50 T0 136

IF termination or Prooram is chosen then end

ELSE IF(SELEDT=999) THEN
END

Hllow user to correct For keyboard input errors

ELSE

MRITECfi); ’PLEHSE SELECT R NUMBER FRON THE MENU’
GO TO 53

ENDIF

Thermal Quantities

NR1TE<¥D ’Input averaQe ambient temperature over

MRITE(#) ’Previous month’

REHDfifi); THVBRM

WRITECWD ’Enter current ambient temPerature’

REHDCt); TCuRHM

NR1TE<¥> ’enter the deep soil temperature For’

NR1TEC$§ ’the area to be evaluated’

REHDCfi), TSDILD

NRITE<#) ’Enter the heat capacity of the air’

REHDfifi), CPRIR

WRITE€#) ’Enter thermal conductivity of the air’

REHDCfi), KHIR

WRITE(#> ’Enter soil-air heat transFer coeFFicient’

REHDfiV), UHIR

NRITE(#) ’Enter thermal conductivits 0? normal soil’

REHD<$>1 KSOILN

NRITEfifi) ’Enter thermal conductivity oF Frozen soil’

REHDfievi KSOILF

NR1TE(*)J ’enter initial moisture content or soil’
RERD(#); IMOIST

NRITEceb, ’enter Final moisture content of soil’

RERDC#), FMOIST

NRITEfit) ’enter heat capacity 0? unfrozen soil

REHDfifibi CPSOLN

NRITEce) ’enter heat capacity of Frozen soil

REHDCfi); CPSOLF

NR1TE¢$D ’enter density of the soil’

REHD(#), DNSOIL

NR1TEC$> ’enter the density or air

REHDcfi), DHHIR

NR1TEK$3 ’enter the thermal diFFusivits o? unfrozen soil
REHDCfl); DSDILH

MRITEfifiD ”enter the thermal diFFusivits of Frozen soil

(:0

m

D

REHDCfibe
NRITEfifl)
REHDCfih,
NRITECfi)
REHDfifi);
HRITEEfi)
fiERDfifi);

Dimensio

NRITEC?)
WRITEifi)
REHDfifib;
NRITECfi)
WRITEiflb
REHDCfi);
NRITEEfi)
REHD(¥);
HRITEfifi)
NRITEC$P
HRITECfi)
MRITECfi)
REHDCfl);
MRITECfi)
HRITE($)
HRITECfi)
REHDC$DJ
WRITEfifl)
WRITEfifi)
HRITECfi)
REHDCfl);
MRITEfifi)
HRITEfifi)
HRITECK)
WRITEifi)
HRITEfifl)
REHDC‘);
WRITEC$3
HRITECfi)
WRITEifi}
MRITEifi)
IF(TNUM

HRITE(*)J

92

DSUILF

’enter tre heat of €usion o? H20
CPFUSE

’enter the heat capacits of water
CPMRTR

’enter the heat caPacite oF ice

CPICE
ns section

’enter the increment size desired in the’
’direction 0? the PiPe (In +eet“’
XINCRE

"enter the number in the

o? increments
I

’direction o? the PiPe (<51)

KHUM

’enter the welocits o? the air desired’
VELCTV

’The time steP is ’KIHCRE/VELCTV’ sec’

’enter the increment size desired in the”

’horizontal direction in the Plane’

’ orthoeonalto the Pipe (In Feet)’

YINCRE

’enter the number

’Horizontal direction in

’to the PiPeCéEi)

VNUM ~

’enter the increment size desired in the’
vertical direction in the Plane '

’orthogonal to the Pipe (In Feet)’

ZINCRE

’enter the number of increments in the

’ vertical direction in the Plane orthogonal’

’to the PiPe ((31: the # oF increments’

OF ifiCPENEfitS ih the
the Plane orthoeonal’

times the increment size should eQual’
’at least 7.5 Feet For accurate results)’

ZNUM

’Enter the increment at which the tube“
’ should be located verticalle’
’increments should extend at least
’Feet below the tube For accurate results’

7 e!
II'J

TNUM
.LT.3 THEN
’Tube increment must be Ereater than 2’

DO TO EBB

ENDIF
WRITE(*)
HRITECfi)
RERD(*)J
HRITECfi)
HI? I TE '3: 1* ,‘J
REHEC‘);

Check InPut

HEITEfifi)
NRITEfifl)
MRITE(#>
HRITEfieb
NR1TE€$D

’Enter the

temPerature diFFerence
’ For stabilite’

TDIFFD

’Enter the maximum number o? time stePsN
to be PerFormed’

TIMER

Quantities

’Rweraoe ambient temPerature’.TRWGHM
’Current ambient temPerature’ TCUPHH

soil temPerature for’,TSDILD
’Thermal conductivity o? the air“ HRIR
’Soil-air heat transfer coefficient’JURIR

’DeeP

MRITE(#)
WRITEfifib
MRITECfl)

NRITEfifi):

NRITEfifi)
MRITEKfl)
NRITEfi!)
NR1TE<$3
MRITEC?)
HRITECfi)
MRITEfifl)
NRITECfi)
NRITECfi)
MRITEC#)
MRITECV)
HRITECfi)
NRITEfiW)
NRITECi)
HRITEE!)
NRITEiib
NEITEfifi)
NRITE(#)
NRITEfifi)
WRITEfifi)
HRITEC¥D
HRITE<$>
MRITE<#}
HRITECfi)
WRITEfififi
HRITECfi)
HRITECW)
HRITEfifi)
NRITE(#)

HRITEC$>
waiTecwa

HRITEEfi)
REHEfifl):

GO TO

J

93

’Thermal conductivity oF normal soil’JKSDILH
’Thermal conductivity of ¢rozen soil’,KSDILF
’Initial moisture content of soil’,IMDIST

’Final moisture content of soil’iFNDIST
’Heat caPacits of the soil’ICPSDLH
’Heat caPacita o? the soil’,CPSDLF
’Densite o? the soil’JDHSDIL
’The increment size desired in the x’
’ direction 0? the PiPe (In Feet)’,XIHCRE
’Thc velocita of the air’,VELCT¥
’The time steP is ’ XINCREHVELCTV ’
’The increment size desired in the’
’horizontal direction in the Plane’
’ orthogonalto the PiPe (In FeetD’JVIHCRE
’The increment size desired in the’
’ vertical direction in the Plane ’
’orthosonal to the PiPe (In feet)’,ZIHCRE
’The number of increments in the
’direction of the eiPe’,KHUM
’The number 0? increments in
’horizontal direction in the
’orthosonal to the PiPe’,VHUM
’The number of increments in the
’wertical direction in the Plane
’to the PiPe’, ZHUM
’The increment at which the tube’

- ~-’
3:?!—

the
Plane

‘0

orthogonal’

’ should be located uerticalls’; TNUM
’The density of air ’JDHRIR

’TemPerature diFFerence resuired ’JTDIFFD
’Maximum number or time stePs ’iTIMEK

of unFrozen soil ’iDSDILN
’.DSDILF

’Thermal di€€usiwit3
’Thermal diffusivity of Frozen soil
’Heat 0? Fusion o? H20 ’JCPFUSE
’Heat capacity of water ’JCPHHTR
’Heat capacity of ice ’JCPICE

’Hre these values correct (V or H)’

HHSNER
IF(HHSNER.HE.¥>
1896

DD TD 16%

94

Rugust 15; 1934 version

Test Data File

mflnlin
G!
0

THVBRM=3B
TCURHM=8
TSDILD=SB
CPHIR=.24
KHIR-.8142
URIR=1.4B
DNRIR=.832?5
KSOILH=S.$
KSOILF=18
IMOIST=.1
FMOIST=.2
CPSOLF=.3
CPSDLH=.2
DNSOIL-95
XINCRE=2
YINCRE=.5
ZINCRE=.5
VELCTV=3688
vuum-za
znum-se
THUM=18
TDIFFD-.81
TINEX=432888
XNUM~25
SRVERH-"H"
PRERNSai
CPFUSE=143.3
CPNHTR=1.8
CPICE=.5
DEBUG=8
DEBUBZ=8
1888 Continue
casesasweawe*eesees##aammummuumwauwwuawmwammsmmumwwmw*mmau
PreProcessor

jdegrees F

;BTUXLBN-F
JETU~FTKHR—FT2~F
3BTU/HR-FT2-F
gLBMXFTS

jFT/HR

flairinflnliflfl

The PrePhocessor establishes a
distribution in the soil based
earth temPeratuhes the average
temperature For the Past month
thermal Properties of the soil.
The actual
superimposed and would include

temperature
on the deep
ambient
and the

distribution is expotential with a sinusoid

additional Factth'

such as short term temperature Fluctuations.
cfififififififlfifififififlflfifififififlfifififimmfifififififlwfiwfififlfififie*#$**#*$$#$*fifi*

IF(FRERHS.EQ.2> THEN
RERD(#> TSDIL

ELSE

Preprocess Rrras

DO 1288 K=BJZNUN+1

TEMPT=TSDILD~C(TSDILD-TRVGRM)¥((FLUHTCK3HFLDHTCZNUM+1)fifiEDD

D0 1188 I=B;RHUM+1
DD 1168 J=~13VHUH+1
TSDILiliJaKaBD=TENPT
CONTINUE

95

1188 CONTINUE
1288 CONTINUE

ENDIF
c
c
c Compute Tube Pcsition
c
TUBEHT=ZNUM-THUH+1
c
c Compute timestep
c
TINSTPzXINCREHVELCTV
5 Compute mass of soil and air elements
- MHSSHR=DHHIR$XINCRE#VINCRE¥ZIHCRE
MHSSOL=DHSOIL$XINCRE$VIHCREfiZIHCRE
c
c Set air temperature distribution in pipe
c

DD 1388 I=8»XHUM+1
TSDIL(I:8;TUBEHT1L3=TCURHM+CCTRVGHM-TCURHM)#
+((FLOHT(I)/FLDRTCHHUM+1)#$2>> .

1888 CONTINUE

2888 Continue

c#*##$####*#*#########*#$#$*#*##fifififlfififlflfiflfifififififififi$§¥§$fi¥fl

c Main Bods

c The heat capacity and moisture content are

c evaluated based on temfierature

caseaweeweaeeae#eeeeeameweeewmeeeeeaeeesemeeeeeeseasnmeeee
L-8

DD 2488 T3 1JTIMEMX
DD 2388 I= IJHNUM
D0 2288 J= BJVNUM
DD 2188 K8 laZNUM

c
c Calculate moisture content as a Function 0? temperature
c Hote moisture content is also a ?unction o? the cooling
c rate and time would be included i? Quantitative
c information were available
c
IF(TSDILCIiJJK,L>.GT.32) THEN
MOICUN=IMDIST+CFMDIST-INDIST)$
+ (1—(ETSUILCI;J,K,L)-32)Hfi(TSDILD-32>$$2)J)
ELSE
MOICDH=FMUIST
ENDIF
c
c Select correct thermal conductivits based on soil temperature
c Select correct thermal capacita based on soil temperature

IF(TSDILEIJJJK,L>.FE.32> THEN
CPSDIL=CPSULH+MUICDH$CPHRTR

KSOILzKBDILH

ELSE IF(TSOIL<I.JJK,L>.LE.38) THEN
CPSUIL=CPSOLF+MDIEOH$CPICE

KSOIL =KSDILF

ELSE
CPSOLCsCPSDLH+MDICDH#CPNHTR+MDICUHfiCPFUSENZ
CPSOILsCPSDLC

fl 1':

nrin

rufifln

filirun

fl

ninrnn

96

KSOIL=KSDILN
END IF

Select appropriate heat transfer coeFFicients

IF element is at the surTace use

IF(K.EQ.ZNUM) THEN
KMRJORflTIMSTP/fiCPSDIL#M SSDL)
K1=KSUIL¥ZINCRE$VINCREHHINCRE
K2=KSDIL¥ZINCRE¥VINCRE/XINCRE
K3=KSUIL#XINCRE*VINCREHZINCRE
K4=URIR*HINCRE¥VINCRE
K5=KSDIL¥ZINCRE$XINCREXVINCRE
K6=KSUIL$ZIHCRE$EINCREHVINCRE

For elements beyond the First two columns

air to soil

but below the surface use soil to soil eQuation

ELSE IF(J.BE.2> THEN

KMHJOR=TIMSTP/(CPSOILfiMHSSDL)
K1=KSDIL#ZINCRE#YINCRE/XINCRE
K2=KSUIL#ZINCREWVINCREXHINCRE
K3=KSDIL*XINCRE*VINCREXZINCRE
K49KSDIL#XINCRE#VINCREXZINCRE
K5=KSOIL§ZINCRE$XINCREHWINCRE
K6=KSDIL§ZINCRE$XINCREfVINCRE

Elements in the First two columns but above or below

the elements surrounding the tube use soil to soil

ELSE IF(K.BE.TUBEHT+2 .UR. K.LE.TUEEHT-2)

KMHJDRBTIMSTP/(DPSDILfiMHSSDL>
K1=KSDIL#ZINCRE*YINCRE/HINCRE
K2=KSDIL$ZINCRE$W1NCREXKINCRE
K3=KSDIL#XINCRE$HINCRE/ZINCRE
K48KSUIL§HINCRE$VINCREKZINCRE
K5=KSDIL#ZINCRE#XINCREHWINCRE
K6=KSDILfiEINCRE#XINCREXVINCRE

Elements adJacent the tube on a dia9onal

ELSE IF(J.E8.1 .HND. K.NE.TUBEHT) THEN

KMHJDR=TIMSTPX(CPSOILfiMHSSDL)
K1=KSDIL§ZINCRE§VINCRE/KINCRE
K2=KSDIL¥ZINCRE$¥INCREXKINCRE
K3=KSDIL#XINCRE#VINCREXZINCRE
K4=KSDIL§XINCRE$VINCEEKZINCRE
K5=KSDIL#ZINCRE¥XINCREXVINCRE
K6=KSDIL$ZINCRE$XINCREEVINCRE

For the element adJacent to the tu
one side air to soil the other sid s

ELSE IF(J.E@.1 .HND. K.EQ.TUBEHT) THEN

KMRJDR=TIMSTP/(CPSOILfiMHSSDL)
K1=KSDIL#ZINCRE$WINCREHKINCRE
K2=fSDIL§ZINCRE$VINCRE/EINCRE
K3=KSDILflXINCRE$VIHCRE/ZINCRE
K4=KSDIL#KINCREfiVINCREXZINCRE

H- 6;

THEN

ntifln

n11

n + ++-+-++»+nonrin

+~++-++!unlfln
6
m
13]

97

K5=UHIR$ZINCRE§HINCRE
K6=KSDIL$ZINCRE#XINCREKVINCRE

For elements directly below the tube use one side air
to soil the other sides soil to soil

ELSE IF(J.EQ.8 .RND. K.E8.TUBEHT-13 THEN
KMHJUR=TIMSTP£CCPSDIL¥MHSSDLD
K1=KSUIL121NCRE*YINCREXHINCRE
K2=KSUIL#ZINCRE#VINCREXKINCRE
K3=KSDIL$XINCREfiVINCREHZINCRE
K4=UHIR$XINCRE$VINCRE
K58KSOIL¥ZINCREfiXINCRE/VINCRE
K6=KSDIL$ZINCRE§HINCREKVINCRE

For elements directly above the tube us” one side air
to soil the other sides soil to soil

ELSE IF(J.EQ.8 .RND. K.EQ.TUBEHT+1) THEN
KMHJDR=TIMSTPHCCPSUILfiMHSSUL)
KlzKSDILfiEINCREfi91NCRE/HINCRE
K2=KSUIL$EINCREfiVINCREXXINCRE
K3=URIR#XINCRE#VINCRE
K4=KSDILWXINCRE#VINCREXZINCRE
K5=KSDIL¥ZINCREfiXINCRE/VINCRE
K6=KSOIL*ZINCREfiXINCRE/VINCRE

element containing tube

ELSE IF(J.ED.8 .RND. K.EQ.TUBEHT)THEN
KMRJOR=TIMSTP/CCPRIRtMRSSHR)
K1=KRIR#ZINCRE#VINCRE/XINCRE
K2=KRIRfiZINCREfiVINCREfXINCRE
K3=UHIR#XINCRE¥VINCRE
K4=UHIR*XINCRE#VINCRE
K5=UHIR*ZINCRE#XINCRE
K6=UHIR#ZINCRE*KINCRE

Convective Heat Transfer EQuation

TSDILCI;J:K;L+1)'
KMHJ0R*(K1#(TSUILCI-11JaKJLD-TSUILCInJ-KJL))+
k2*':T'SIjIL-"1+1I'JIksJLD~TSIJIL<II~TJK)L—) }+
K3$€T301L£IJJJK~1JLD~TSDILCIiJaKiL35+
K4$<TSDILCIiJsK+1;L)~TSDILCIsJikaLDJ+
K5*(TSDILCIaJ-19KsL)-TSUIL(IiJiK1L)3+
K6¥CTSUIL(I;J+1:K;L)~TSUILCI;J,K,L)D3+
T'SOILCI-l I J/KJL.)

GO TO 2188

ENDIF

H .t transFer eQuation

cl)
0

SDILKIJJJKJL+138
MRJGRflfiKififiTSDILfiI-i,JJKJL)~TSDIL(IJJJKJL3>+
K2$iTSOILC1+1aJiK2LD~TSDIL(I,J,K,L)3+
KafifiTSUILCI,J,K-1»L)-TSDIL<IJJiKJL))+
K4§fiTSDIL<IaJiK+1JLb~TSDILci,JJKJL))+
K5*(TSDILKIJJ-19K3L3~TSDILEI,J,K;L))+

T
K

fltfironrlfl +-+

FIG

ntnrin

ntirun

flflaifln

nlfiflfl nrin

nrnn

98

K6¥QTSOILCIJJ+1aKiL)-TSUILCI;JJKJLDb)+
TSOILCIaJiKiLfi

HdJust For temperatures crossing state change
boundries

IF the element is the air Pipe no change is necessara

IF(J.EQ.8 .HND. K.ED.TUBEHT) 88 TD 2188

.~ _-

I? the initial temperature is above ea
IF(TSOILfiliJiKJL).BE.32) THEN

find the ending temperatures is above 32
degrees then no adJustment is necessars

IF(TSDILcIJJIKJL+1)QUE-32) THEN
80 TD 2188

find the Final temperature is below 38 then a comPound
adjustment is necessars

ELSE IF(TSDILfiI,J,KJL+1>.LE.38> THEN
DELTHiaTSDIL<IiJiKiLD~TSUILCIiJJKiL+13
Q-CPSOLNfinELTHI

eiac SDLN#(TSUIL<I,J,K,L~32)
o2=cpsoLcs2

IF(Q-Qi .LT. D2) THEN
DELT8228-81HCPSDLC
TSDILcI,JiK,L+1)=32~DELTH2

ELSE

DELTH3=<Q-Qi-Q2)/CPSULN
TSOILcIiJiKiL+1>=38-DELTHS

I? the Final temperature is below 32 but above 38 the
temperature reduction is reduced to account For the
higher soil heat capacity below 32

ELSE IF(TSUILfiliJJK;L+1).LT.32 .RND .PT. 38} THEN
DELTHI=32~TSDILCI:JsK;L+1)
DELTRZ=CPSDLN/CPSDLC#DELTR1
TSDIL£IiJiKiL+1)=S2-DELTHE

ENDIF

IF the beginning temperature is below 38 degrees
ELSE IF(TSDILCIJJJKJL).LE.38) THEN

find the ending temperature is below 38
degrees then no adJustment is necessars

IF(TSDILEIIJJKJL+1}.LE.38) THEN
58 TD 2188

. compound

(d
N
‘+
w-
.J
n
J
.J
(u

find the Final temperature is above ”
adJustment must be made

nrnn

niin

nt1r:n:1

fl 0 fl fl

n11

nrun

+

99

ELSE IF(TSOILCIaJiK,L+1).OE.32> THEN
DELTRI=TSOIL<I,J,K,L+1D—TBOILCIJJJK,L)
OsCPSOLFenELTHI
O1=CPSOLF#<38~TSOIL<IaJiKiLD
OQeCPSOLC$2

IF(O-Oi. LT. O2) THEN

DELTRE =o-O1'cps OLI:

TSOILCI J K L+13=38+DELTH2

ELSE

DELTR3=£O~O1~O£>£OPSOLN
TSOILKIJJJKiL+1>=32+DELTHS

find iF the Final temperature is between 38 and 32 degrees
then the Final temperature is reduced

ELSE

DELTRI=TSDIL<IJJ K. L+1)-E '8
DELTR2 =3P'SDLF.ILPS ULCEfiDELTHi
TSDIL(1.MT K.L+1)=38+DELTH2
ENDIF

IF the beginning temperature is between 38 and 32

ELSE

and the Final temperature is above 32 degrees the Final
temperature is increased

IF(T 381L£I J. K L+1>. BE. 32) THEN

DELTRi-TSOILCI..J K.L+1)-‘2

DELTR2=CPSOLC(:PEDLNtDELTHl

TSUILfiiiJiiiL+iD=TSDILfiI J. K L+1>+DELTHZ

and the Final temperature is below 38 degrees then the

temperature is Further reduced

ELSE IF(TSUIL(IaJ:KJL+1).LE.S8) THEN
DELT81="8- TSDILC IJJJKJL+1)
DELTH2=CPS OLC '8 'SDLFfiflELTHI
T'SDIL(I; J.K;L+1)=S8-DELTH2

find the Final temperature is between 38 and 32
then no adjustment is necessary

ELSE

GO TO 2188

ENDIF
ENDIF
CONTINUE
CONTINUE
CONTINUE

Test For temperature constant condition
IF temperature is less than the temperature chan
desired Jump to program Finish

q;-
II!

WRITE(#)J TSDILfikNUN 8.- TUBEHT L+1). ’Jews we ’
IF(TSDILC.4.NUM. 8 TUBEHT L+1J-T-UIL'VP|1 8 TUI’EHT. L)
LT. TDIFFD) 88 TD 2? 88

nlnrnnri

ntwnim
U) ‘4!
Ltd
6C9

fiflfl

FIfiiTth

mrhnru

(on

' fl +-++

can

100

Taking advantage of the plane of semmetra around the
tube in the x-: plane Fill TSDILfiIJ~1JKiL+1> arras
with TSDILCIJ+IaKiL+1) values

DO 2328 M=1JNNUH

D0 2318 N=1aZNUN
TSOILfifli-l:NJL+1)=TSOIL<M;laNsL+1)
CONTINUE

CONTINUE

Match end wall temperatures to array end temperature
in order to satisfy insulation condition

DO 2348 J81;VNUN

DO 2338 K=IJZNUN
TSOILCBIJ;KaL+1)=TSOIL(13J»K:L+1)
TSDILfiXHUM+IJJJK3L+1D=TSDIL<XHUMJJJKJL+1>
CONTINUE

CONTINUE

Reset tube inlet temperature to current ambient
TSDILIBIBJTUBEHT,L+1)=TCURRM
Move L+1 temperatures into the L arras For next loop

DO 2398 D=8;NNUM+1

D0 2398 B=~IJVNUM

DO 2398 C=112NUN
TSOILCD:BICJLD=TSOIL(DJBJC;L+1)
CONTINUE

CONTINUE

IF maximum time step is reached print out message

WRITEC*); ’Maximum time exceeded’
80 TO 2888

IF stable temperature is achieved print message

HRITEfifi); ’Output temperature stable at’
MRITEfifi). TSOILiNNUMJBJTUBEHT.L+1); ’de9rees.’

Hsk i? user wishes to save parameters

MRITEfiflD, ’Do sou wish to save the parameters?’ (9 or N)
REHECfi); PRRHN

IFQPHRHN.EQ.’V’)THEN
NRITEfifiifisERR=4488§THVORM,TCURHNJTSDILD:CPHIRJ
KHIRJUHlfiinNRIR,~SOILNaKSOILFaIMOISTJFMOISTitPSOLF;
CPSOLN,DNSOILJNINCREJVINCRE,ZINCREJVELCTV,VNUM,ZNUM’
TNUM;TDIFFD,TIMEX,XNUH;$HVEHN.PRERNSJCPFUSEJCFNRTRTOPIOE
ENDIF

Rsk i? user wishes to store the array walu

I7!

3.

NRITE(#3, ’Do sou wish to save th
REflnfifi), SRVERN
IF(SRVERN.EO.’9’>THEN

temperatures (V or N)

11;

(.0 is)
(5' SI
‘5' G!

4188

4488

4588

4688

101

NRITECSI¥JERR=4S88D TSOIL
ENDIF

CLOSECoxERR=4688D
CLOSEfiéiERR=4788D

NRITEfifi); ’NORNHL TERMINHTION’
STOP

Error messases and recovers

NRITE<#) ’Open Failed For parameter File’

MRITECfiD ’Do sou wish to attempt to solve the problem’
HRITE(#> ’ and trs asain?’

NR1TE€$> ’Oorrect error cause and answer 9’

MRITE<$> ’IF sou wish to Quit enter N’

RERD(#> H

IF(R.EO.’V’) GO TO 18

STOP

HRITECfi) ’Open Failed For data File’

MRITECSJ ’Do sou wish to attempt to solve the problem’
HRITECfi) ’ and trs asain?’

NRITE£$D ’Correct error cause and answer s
WRITE(#> ’IP sou wish to Quit enter N’
REHDCfi) R

IF(H.EO.CV’) GO TO 28

STOP

HRITE(*) ’Read Failed For parameter File’
NRITEfifi) ’Do sou wish to attempt to solve the problem’
HRITE(*) ’ and trs asain?’

MRITECfi) ’Correct error cause and answer 9’

WRITEESD ’IF sou wish to Quit enter N’

REHDCfi) H

IF(R.EO.’V’> GO TO 78

STOP

NRITE<#) ’Read Failed For data File’

NPITECfi) ’Do sou wish to attempt to solve the problem’
NRITECS) ’ and trs again?’

HRITEC#) ’Oorrect error cause and answer V’

NR1TE£¥> ’IF sou wish to Quit enter N’

REHD(#) H

IF(H.EO.’V’) GO TO 1188

STOP

NRITEfia) ’Nrite Failed For parameter File’

MRITECS) ’Do sou wish to attempt to solve the problem’
MRITECfi) ’ and trs again?’

MRITEC$> ’Correct error cause and answer V’

NRITEC#D ’IF sou wish to Quit enter N’

RERDC#) H

IF(H.EO.’V’) GO TO 2758

STOP

MRITECS) ’Nrite Failed For data Pi
NRITE<#) ’Do sou wish to attempt t
NR1TE<$D ’ and trs again?’
NRITEifib ’Correct error cause and answer 9’
NR1TE£§> ’IF sou wish to suit enter N’
PEHDCSD R

IF(R.EO.’V’) GO TO 2388

STOP

NRITECfi) ’Olose Failed For paramet r 1
HEITEflfi) ’Do sou wish to attempt to sole
NRITEifib ’ and trs asain?’

' I

I

le'
o solve the problem’

57}

the problem’

102

[ll

NRITEfifib ’Correct error cause and answer 7
NRITEifi) ’IF sou wish to Quit enter N’

REHDfifl) H
IF(H.EQ.’V’) GO TO 3488
STOP

NRITE(#> ’Close Failed For data File’

MRITE(*) ’Do sou wish to attempt to solve the problem’
NR1TEC¥> ’ and trs asain?’

NRITEce) ’Correct error cause and answer s!
NRITEcwb ’IF sou wish to suit enter N’
RERD(#> H

IF(H.EO.’V’) GO TO 3688

STOP

END

BIBLIOGRAPHY

10)

11)

12)

13)

BIBLIOGRAPHY

Heap, R. D. Heat Pumps. A Halstead Press Book, Copyright
1979 by R. D. Heap, page 11.

 

McGuigan, Dermot, and McGuigan, Amanda. Heat Pumps, An Efficient
Heating & Coolinnglternative. Charlotte, Vermont: Garden Nay

 

 

Publishing, 1981, page 29.
Reference 2, page 53.
Reference 1, page 25.

Kirschbaum, H. S. and Veyo, S. E. "An Investigation of Methods
to Improve Heat Pump Performance and Reliability in a Northern
Climate." EPRI EM-319 Project 544-1 Final Report, Prepared by
Westinghouse Electric Corporation, Pittsburgh, Pennsylvania,
January 1977.

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.
OakTRidge, Tennessee: Oak Ridge National Laboratory, National
Technical Information Service, 1981, Chapter 4.

Gannon, Rogert. "Ground-water Heat Pumps." Popular Science, 1978.

 

Bylinsky, Gene. "Water to Burn." Fortune, October 20, 1980.

Sumner, John A. Domestic Heat Pumpg. Chalmington, Dorchester:
Prism Press, Stabel Court. Dorset DT2 OHB, Great Britain.
Copyright 1976 by John Sumner, page 36.

 

Metz, Phillip D. "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.

Bose, J. E. "Earth Coil/Heat Pump Research at Oklahoma State
University." Stillwater, Oklahoma: School of Technology,
Oklahoma State University.

Nasserman, M. "Crawl-Space Assisted Heat Pump with Numerical
Model.“ M.S. thesis, University of Tennessee, 1983.

Robb, Ronald J. "An Analytic Analysis of the Air-to-Air Ground
Coupled Heat Pump." March 9, 1983.

MN
4500

104

Reference 11.

Reference 9, page 41.

Reference 11.

Jennings, Burgess H. Environmental Engineering,_Analysis and

Practice. New York, New York: International Textbook Company,
I970, page 142.

 

Reference 17, page 139.

Koorevaar, P., Menelik, G., and Dirksen, C. Elements of Soil
Physics. New York, New York: Elsevier Science Phblishers, 1983,
Equation 4.31.

Tsytovich, N. A. The Mechanics of Frozen Ground. McGraw-Hill
Book Company, 1975, pageg34.

 

Yong, Raymond N., and Warkentin, Benno P. Soil Properties and
Behaviour. New York, New York: Elsevier Scientific Publishing
Company, 1975, page 417.

 

Reference 20, page 70.
Reference 20, Section 11.11.
Reference 7.

Boyer, Richard, and Savageau, David. Places Rated Almanac.
New York, New York: Rand McNally, 1981, page 26.

Reference 17, page 136.
Reference 17, page 164.

Potter, Merle C., and Foss, John F. Fluid Mechanics. New York,
New York: The Ronald Press Company, 1975.

 

Reference 17, page 138.
Reference 17, Table 4.3.
Manufacturers data.
Reference 3, page 41.
Consumers Power, 1984 rates.

Reference 3, page 23.

105

35) Aero Rental, 1984 rates.
36) Reference 17, page 422.
37) Grainger Fan and Blower Catalog, Fall 1983 edition.

38) Jumikis, Alfred R. Thermal Geotechnics. New Brunswick, New
Jersey: Rutgers University Press, 1977.

 

39) Reference 11.
40) Reference 11.

General References

Potter, Merle C. Numerical Methods in the Physical Sciences. Englewood
Cliffs, New Jersey: Prentice—Hall, Inc., 1978.

 

Holman, J. P. Heat Transfer. New York, New York: McGraw-Hill, Inc.,
1976.

 

Jumikis, Alfred R. Soil Mechanics. Princeton, New Jersey: 0. Van
Nostrand Company, Inc., 19621

 

Jumikis, Alfred R. Thermal Soil Mechanics. New Brunswick, New Jersey:
Rutgers University Press, 1966.

 

Reedy, w. R., and Bullock, C. E. ”Northern Climate Heat Pump Field
Performance Evaluation." EPRI EM-2319. Carrier Corporation,
July 1982.

Reay, D. A., and Macmichael, D. B. A. Heat Pumps Design and Application.
New York, New York: Pergamon Press, 1979, page 39.

Heins, Royal 0., and Rotz, C. Alan. "Final Report to Detroit Edison
on Feasibility of Greenhouse Production in the Detroit Area
Utilizing Heat Pumps for Temperature Control." September 1981.

"Heat Pump Demand Characteristics, A Study of the Impacts of Single-
Family Residential Heat Pump Technologies on Electric Utility
System Loads." EPRI EA-2074, Project 1100-1, Final Report,

October 1981, Prepared by Gordian Associates, Inc., Washington, DC.

Collie, M. J., ed. Heat Pump Technology for Saving Energy. Park Ridge,
New Jersey: Noyes Data Corporation, 1979.

 

"Backyards Heat New York Homes." The Lansing State Journal, February 20,
1983.

Faltermayer, Edmund. "Solar Energy Goes Underground.” Fortune.
December 27, 1982.

106

Andrews, J. N., "Ground Coupled Solar Heat Pumps: Analysis of Four
Options." ASME Solar Energy Division Conference, Reno, Nevada,
April 27-May 1, 1981, Department of Energy and Environment,
Brookhaven National Laboratory, Upton, New York.

 

 

M11111
3

  

   
 

   

TAYE U

S

1293 0

 

 

 

 

 

 

 

 

 

 

 

           

1111

 

   

      

   

 

TYLB

11111“

716