EXPERIMENTAL AND COMPUTATIONAL INVESTIGATIONı
OF SNOW MELTING ON A HYDRONICALLYı
HEATED CONCRETE SLABı
Byı
DIEGO PATRICIO ESPIN TOBARı
Bachelor of Scienceı
Anny Polytechnic Schoolı
Quito, Ecuadorı
1998ı
Submitted to the Faculty oftheı
Graduate College of theı
Oklahoma State Universityı
in partial :fulfillment ofı
the requirements forı
the Degree ofı
MASTER OF SCIENCEı
December, 2003ı
EXPERIMENTAL AND COMPUTATIONAL INVESTIGATIONı
OF SNOW MELTING ON A HYDRONICALLYı
HEATED CONCRETE SLABı
Thesis Approved:
ii
ACKNOWLEDGEMENTS
This thesis could not been completed without the participation ofse;veral persons.
I wish to express my gratitude to all those who directly and indirectly participated in this
process.
In addition, I would like to thank my advisor Dr. Spitler for his practical and
theoretical guidance and support. Dr. Spitler has used his expertise in the field of thennal
systems to guide me in the development ofthis thesis.
The love and support from my family have been ofmajor importance. I would
like to thank my wife" Siri, for all her support and understanding. She has always been
there for me. Further I would like to thank my "adopted" American grandparents for aU
their support and love during my stay in Stillwater. Last but not least, to my parents
(Lupita and Ruben Espfn) for their great support throughout my life. This has enabled me
to be where I am today.
My gratitude to Dr. Mark Rockley who has been very supportive with his critical
questioning and for the long hours spent helping with the revision ofthis thesis.
I would further like to thank to my fellow research assistants and friends.
Xiaobing Liu and Dongyi Xiao for their advice and help during the development ofmy
research.
iii
r---__- ---------------------/'
TABLE OF CONTENTS
CHA..PTER 1 INTRODUCTION _ 1.1 OVERVIEW 1.2 STEADY STAIT AND TRANSIENT MODELS 1.3 PROBLEM STATEMENT 1.4 OBJECTIVE J CHAPTER 2 - LITERATURE REVIEW 2.1 INTRODUCTION 2.2 HYDRONIC SNOW MELTING SYSTEMS OVERVIEW 2.3 STEADY STATE MODELS 2.3.1 1-D Steady State Models 2.3.2 2-D steady state models 2.4 TRANSlENT MODELS '" 2.4.1 I-D transient models 2.4.2 2-D transient models 2.5 EXPERIMENTAL INVESTIGATIONS 2.5.1 Snow melting by heating from the bottom (Aoki, et al. 1987) 2.5.2 Hockersmith experimental investigation CHA..PTER 3 - -..........•...................................... IV
CH.APTER ' 'PAGE
EXPERIMENTAL APPARATUS..............•...........•.............•................•....................... 3.1 INTRODUCTION 3.2 CONCRETE SLAB 3.2.1 Base 3.2.2 Heating Systen1 3.2.3 Instrumentation 3.3 SNOW MAKING / ENVIliONMENTAL CHAMBER 3.3.1 Top section 3.3.2 Central section 3.3.3 Bottom section 3.4 MECHANICAL COOLING 3.5 SNOW MAKING EQUIPMENT 3.6 DATA ACQUISITION CONTROLLER 3.7 CONCRETE THERMAL CONDUCTION COEFFICIENT MEASUREMENT
APPARATUS 3.8 SPECIFIC HEAT CALCULATION APPARATUS CHA.PTER 4 EXPERIMENTATION AND RESULTS 4.1 INTRODUCTION 4.2 THERMAL PROPERTY MEASUREMENTS 4.2.1 Density measuren1ent v
CHAPTER PAGE
4.2.2 Specific heat of concrete
4.2.3 Thermal conductivity
4.2.4 Water diffusion on concrete experiment..
4.3 MATHEMATICAL MODEL OF THE CONCRETE SLAB 4.3.1 Heat conduction 4.3.2 Specific heat and density 4.4 INITIAL CONDITIONS AND PHYSICAL PROPERTIES AS REQUIRED INPUT 4.4.1 Experimental procedure 4.4.2 Power input calculation 4.4.3 Initial conditions 4.5 EXPERI1VfENTS 4.5.1 Dry surface condition 4.5.2 Snow-melting experiments 4.5.3 Revised model with moisture penetration 4.6 SUMMARy CHAPT'ER 5 CONCLUSIONS AND RECOMENDATIONS 5.1 CONCLUSIONS 5.2 REcoMMENDATIONS
BffiLIOGRA.-IT ,
APPENDIX A :
-~-.---:---:--::-:::-::--::------~---------------------
CHAPTER . PAGE '
INSTRUMENT CALmRATION •.........•••....•.•••. ,
A.l INTRODUCTION A.2 VARIABLES & CALIBRATION , ." APPENDIX B DATA LOGGER CONTROLLER B.l INTRODUCTION B.2 DESCRIPTION OF THE SYSTEM APPENDIX C GEOMETRIC INFORMATION OF CONCRETE SLAB C.1 CONCRETE SLAB C.2 MAIN THERMOCOUPLE PLACEMENT DIMENSIONS C.3 CENTRAL SECTION APPENDIX D ...•.........•.......•..•..•••......•....•... .•....•......•..••••••••.•.•..•••..•••..•.••.••••..•.•.• •... CHARTS AND DATA USED FOR CONCRETE CONDUCTIVITY
CALCULATION ; D.I SAMPLE 1: D.2 DATA OBTAINED FOR SAMPLENUMBER2: vu
TABLE OF FIGURES
FIGURE 2.2-1 : SERPENTINE PIPE CONFIGURATION IN A HYDRON1CALLY- HEATED SLAB IN
PLANE AND CROSS SECTIONAL VIEW........ ...............•.......••• " FIGURE 2.4-3: MODEL DOMAIN SHOWING FINITE DIFFERENCE GRID AND BOUNDARY
FIGURE 2.4-4: THE LOCAL COORDINATE SYSTEM ATTHE EASTFACE OFA TYPICAL FINITE
FIGURE 2.4-8: SCHEMATIC REPRESENTATION OF HEAT rRANSFER IN THE THREE-NODE
FIGURE 3.2-1: THERMOCOUPLE AND TUBING PLACEMENT DURING CONCRETE POURING
FIGURE 3.2-6: A) AUXILIARY EQUIPMENT. Bj AUXILIARYEQUIPMENT COVERED WITH
FIGURE 2.2-2: SLINKY PIPE CONFIGURATION IN A HYDRONICALLY-HEATED SLAB ., CONDITIONS VOLUME CELL : " FIGURE 2.4-5: BLOCK DEFINITIONS ON A TYPICAL GRID LA YOUT FIGURE 2.4-6: GRID GENERATED FOR HYDRONIC SYSTEM .' FIGURE 2.4-7: MASS TRANSFER ON SNOW MELTING MODEL. " SNOWMELTMODEL PROCESS FIGURE 3.2-2: CONSTRUCTION OF CONCRETE SLAB FIGURE 3.2-3: CONCRETE SAMPLES DURING SLAB CONSTRUCTION FIGURE 3.2-4: CONCRETE SLAB AND BASE FIGURE 3.2-5: CONCRETE SLAB SIMPLIFIED PIPING DIAGRAM INSULATION INSIDE THE SNOW CHAMBER FIGURE 3.2-7: POSITION OF THERMOCOUPLE AT THE CONCRETE SLAB FIGURE 3.3-8: ENVIRONMENTAL CHAMBER SECTIONS viii
FIGURE PAGE
FIGURE 3.3-9: A) LENS OPENING. B) DIGITAL CAMERA SPACE IN ALTERNATE PLUG FIGURE 3.5-12: A) MODIFIED NOZZLE SCHEMATICS, B) NOZZLES OUTSIDE CHAMBER, C)
FIGURE 4.3-3: CENTRAL SECTION. A) CONSTRUCTION MODEL, B) MA THEMATICAL
FIGURE 4.5-7: AVERAGE SURFACE TEMPERATURES WITH DlFFERENT VALUES OF CpFOR
FIGURE 3.3-10: MECHANICAL COOLING DIFFUSER FIGURE 3.4-11: SNOWMAKING MACHINE SCHEMATICS NOZZLES INSIDE THE CHAMBER FIGURE 3.7-13: CONCRETE TESTING PACK FIGURE 3.8-14: CONCRETE SAMPLE FOR CP EXPERlMENTATION FIGURE 3.8-15: CALORIMETER INSULATION FIGURE 3.8-16: DIGITAL BALANCE FIGURE 4.2-1 : SCHEMATIC OFHEAT CONDUCTIONON CONCRETE SAMPLE FIGURE 4.2-2: WATER DIFFUSIONIN CONCRETE SLAB SAMPLE REPRESENTATION FIGURE 4.3-4: RESISTANCE REPRESENTATION OF SLAB'S BASE FIGURE 4.5-5: A VERAGE TEMPERATURE COMPARISON IN DRY CASE SCENARIO FIGURE 4.5-6: AVERAGE SURFACE TEMPERATURE FOR SNOW EXPERIMENT CONCRETE FIGURE 4.5-8: SNOWFREE AREA RATIO VS. TIME FOR SM] EXPERIMENT FIGURE 4.5-9: FD-RG PREDICTED SNOWFREEAREA RATID 0.00, MEASURED: 0.00 FIGURE 4.5-1 0: FD-RG PREDICTED SNOW FREE AREA RATIO 0.08, MEASURED: 0.02 FIGURE 4.5-11: FD-RG PREDICTED SNOWFREE AREA RATIO 0.25, MEASURED: 0.20 ix
FIGURE PAGE
FIGURE 4.5-12: FD-RG PREDICTED SNOWFREE AREA RATIO: 0.54, MEASURED: 0.48 .•..... FIGURE 4.5-13: FD-RG PREDICTED SNOW SURFACE FREE AREA RATIO: 0.71,
FIGURE 4.5-14: FD-RG PREDICTED SNOW SURFACE FREE AREA RATIO: 1.00,
FIGURE 4.5-15: FD-RG PREDICTED SNOWSURFACE FREE AREA RATIO: 1.00,
FIGURE 4.5-23: SURFACE TEMPERATURE COMPARiSONBETWEEN MODIFIED AND
FIGURE 4.5-24: SURFACE TEMPERATURE COMPARISONBETWEEN MODIFIED AND
FIGURE A-I: CALIBRA TION CHART FOR THERMOCOUPLE INSIDE SAMPLE FOR CP
FIGURE A-2: CALIBRATION CHARTFOR THERMOCOUPLE OUTSIDE OF SAMPLE FOR C p
MEASURED: 0.65 MEASURED: 0.95 MEASURED: 1. 00 FIGURE 4.5-16: FD-RG MODEL AND FV-BFG MODEL COMPARiSON FIGURE 4.5-17: SNOWFREE AREA RATIO VS. TIME FOR SM2 EXPERIMENT • ........................... FIGURE 4.5-18: FD-RG PREDICTED SNOWFREE AREA RATIO 0.00, MEASURED: 0.00 FIGURE 4.5-19: FD-RG PREDICTED SNOWFREE AREA RATIO 1.00, MEASURED;' 0.95 FIGURE 4.5-20: FD-RG PREDICTED SNOW FREE AREA RA TIO 1.00, MEASURED: 1.00 FIGURE 4.5-21: ESTIMATED MOISTURE PENETRA TION DEPTHIN SM1 EXPERIMENT FIGURE 4.5-22: ESTIMATED MOISTURE PENETRATION DEPTH IN8M2 EXPERIMENT ORIGINAL FD-RG MODEL INEXPERIMENTSM1 ORIGINAL FD-RG MODEL IN EXPERIMENT8M2 DETERMINATION DETERMINATION x
FIGURE , PAGE
FIGURE A-3: CALIBRATION CHARTFOR THERMOCOUPLE IN WATER BATHFOR Cp
DETERMINATION _ , 155
FIGURE A-4: CALIBRATION CHARTFOR THERMOCOUPLE INTO CALORIMETER FOR Cp
FIGURE 0-2: POWER MEASUREMENT ON THERMAL CONDUCTION CALCULATIONFOR
, FIGURE 0-4: POWER MEASUREMENTFOR THERMAL CONDUCTION CALCULATION FOR
DETERMINATION 155
FIGURE A-5: CALiBRATION CHARTFOR THERMOCOUPLE ATINLET OFHYDRONIC PIPING 156
FIGURE A-6: CALIBRATION CHARTFOR THERMOCOUPLE ATOUTLET OF HYDRONIC PIPING 156
FIGURE B-7: TIME BASE DIAGRAM 158
FIGURE B-8: LIGHTCONTROL CIRCUIT 160
FIGURE 8-9: LISTOFSOURCE CODE FOR BASIC STAMP 11.. 161
FIGURE B-1O: COMPLETE DIAGRAMOFDATA ACQUISiTiONCONTROLLER SYSTEM 162
FIGURE C-I: UPPER SLAB SURFACE 163
FIGURE C-2: SLAB FRONTAL DIAGRAM 164
FIGURE C-3: CENTRAL SECTION OF SLAB 164
FIGURE C-4: MAIN THERMOCOUPLE PLACEMENT 165
FIGURE C-5: CENTRAL uREAL" AND "MODELED" SECTION DIMENSIONING 166
FIGURE D-1: TEMPERA TURE CHART IN SAMPLE 1 _ 167
SAMPLE 1 168
FIGURE 0-3: TEMPERATURE CHARTIN SAMPLE ·· 169
SAMPLE 2 167
xi
CHAPTER]
INTRODUCTION
1.1 Overview
Freezing rain, frost, and snow formed by water or ice crystals falling and
collecting on the ground can create dangerous travel situations. Snow, freezing rain or
frost on airport runways, on highways and roads, or on bridges can be a serious threat to
human travel safety. Bridges are especially critical to travel safety - not being in contact
with the ground, they freeze quicker than the adjacent roadway.
To increase travel safety, snow can be removed or melted from paved surfaces.
The use of chemicals (e.g. salt) or mechanical devices (e.g. snow plows) is a frequent
practice. The most commonly used chemical on roads is salt (Sodium Chloride). Salt
interacts with water and ice and decreases the freezing point. Salt is an inexpensive
deicing agent and it is efficient in melting snow and ice except When the surface
temperature is lower than -9°C (1 5°F).
The primary disadvantage ofusing salt is its corrosive interaction with the bridge
deck reinforcing steel and bridge structural steel over the long term. The resulting repair
costs ofbridges are often very high. The salt will further interact witb the vehicles and
over time corrode them also. Other disadvantages of chemical Use are: environmental
.J
I
damage on Vegetation and adjacent water streamS (LOfgren 2001) and the difficulty of
application before a snowstorm arrives. Salt is not applied before storms because it will
not react without moisture and can be swept away before the snow anives.
Another method of ice/snow removal is the use of mechanical devices (e.g. snow
plows), but these can damage the pavement and can be associated with high maintenance
costs. Machinery vibration and improper use oftbe mechanical device lead to surface
and structural damage (Nixon 1993). Again, it is a purely reactive method is used once
the freezing condition is present and t1;terefore the safety concern is present. Moreover,
severe environmental conditions may preclude use ofthe equipment. For these reasons,
alternative ice melting methods that permit cost savings, lower maintenance, and reduce
road and bridge damage should be considered.
One alternative method is to embed hydronic tubing in critical points ofroads,
highways and bridges to melt the snow. A fluid can be heated using a ground source heat
pump that transfers heat from the ground to the fluid. The warm fluid is circulated
through the embedded pipes in order to prevent freezing conditions on the pavement. In
this way, the pavement surface can be kept free from ice or snow. This system is being
investigated by researchers at Oklahoma State University and is referred to as the "Smart
Bridge Project". The name comes from the use of local and remote weather data to
forecast icing conditions and automatically control the heating system. The objective is maintain a surface free ofice or snow during the whole storm period.
The key problem for the Smart Bridge System is the high initial cost. An
important eIl1phasis of the research at OSU is the optimal design and control ofthe Smart
2
Bridge system. Reliable mathematical models of all the system components are required
in order to accomplish this task. The OSU research group has developed new component
models ofheat pumps, ground loop heat exchangers, and bridge decks for use in the
modular simulation environments, TRNSYS and HVACSIM+. Testing and validation the component models are essential to ensure their accuracy and the validity ofthe
system simulation. This thesis is concerned with validation ofthe bridge deck model.
1.2 Steady state and transient models
Snow melting models can be classified as either steady state or transient. Earlier
analyses of snow-melting systems perfOlmance and guidelines for the design of such
systems have been based upon steady-state oonditions.'In a steady state analysis, the
response ofthe dynamically heated slab and rapidly changing moisture conditions are considered (i.e. no considerations are taken from the storm history). Since storms are
dynamic in nature, steady state models can only give a very TOUgh idea of the actual
system behavior. Moreover, the system capacity can be over (or under) estimated and
therefore make the system inadequate for the application, or unnecessarily expensive.
By contrast, transient models take into account the dynamic response of the slab
to changing weather conditions. Two transient models h.ave been developed to analyze
snow melting. Chiasson, et at (1999) developed the first model used as part of the "Smart
Bridge" project at oklahoma State University. Subsequently, Liu, et al. (2002) improved
upon this early model by refining the correlations and grid corresponding to the pipe
approximation used by the previous model. Rees, et al. (2002) developed an alternative
model that accounted for the mixed geometries on the grid model used. This model uses 3
different mathematical approach to solve the heat balance due to the transient conditions.
However, the researchers have also established the need of further refinement and
experimental validation for these models.
1.3 Problem statement
An important requirement in the design of a snow melting system is the
determination ofthe optimum design parameters for a heated slab. These design
parameters which must be optimized are: the depth and spacing of the heating elements
within the slab, the operating temperatUre ofthe working fluid, and the proper layout of
the embedded piping in order to obtain the desired surface condition (i.e. ice free). These
parameters can only be optimized with a reliable model capable ofpredicting the
response of the system to highly transient storm conditions.
Both ofthe models actually used for simulation require experimental validation a controned environment. Experimental (field) validation has been performed as part of
the "Smart Bridge" project. However, the environmental conditions could not be
controlled, and it is desirable to perform experiments under well-controlled conditions,
where long and short wave radiation losses can be quantified more accurately. The results
of the experimental validation used in the "Smart Bridge" project showed a difference of
up to 2 °C (3.6 OF) from the measured temperature on bridge deck surface during the
testing period (Liu, et al. 2002). The model presented in ASHRAE I090-RP was partially
validated by Hockersmith (2002). He only tested the snow-melting portion ofthe
algorithm. Moreover in his validation he used an electric heated plate that instantly
transferred the power to the surface. The actual transient response of the main component
4
(the slab) on a bridge deck was not considered. A real concrete slab needs to be tested
under controlled conditions.
1.4 Objective
The objectives ofthis study are to experimentally validate the two transient, twodimensional
models of a hydronically-heated concrete slab with snow melting boundary
conditions, under an environmentally controlled condition.
5
CHAPTER 2
LITERATURE REVIEW
2.1 Introduction
In this section, an overview ofJ;1ydronic snow melting system technology is
presented! first. This overview is followed by a review of previous published research
related to modeling of snow-melting systems. The revie,wed models are presented in two
groups. The first group ofmodels uses steady state calculations. In the second group,
models based on transient analysis tools are reviewed. Finally, experirnental
investigations are discussed.
2.2 Bydronic snow melting systems overview
Practical snow melting systems use embedded piping or electric heating elements.
Both types of systems maybe used to reduce safety concerns over bridges and roads and
to keep these surfaces maintenance free in the winter. The maintenance usually
associated with winter periods includes plowing and applying de-icing chemicals. The
use of snow melting systems can eliminate all this work and prevent potential damage the concrete and structures by snow removal equipment and chemicals. Furthermore, it
can prevent snow accumulation.
6
~---------------------~~-------
heating element. In the case ofbydronic systems, an embedded pipe circulating a mixture
of antifreeze (e.g. propylene glycol or ethylene glycol) and water is used as the heating
element. In both cases "the system time constant is in the order of hours" (Rees. et al .•
2002). This is caused by the large thermal mass of the system and is important because the highly transient nature of weather.
There are two configurations ofhydronically-heated pavement systems that have
been reported. The first one consists of pipes embedded in the pavement material evenly
spaced and connected with V-Shaped tubes as shown in Figure 2.2-1. This system is
commonly used for snow melting applications.
,,-...., ,,-... /""0. ,,-...., ,,-... ,.,--..
, I
'--"
I 000 0 0 000 0 000
Figure 2.2-1: Serpelltille pipe configuration in a hydronically- heated slab in plane and cross sectional
view
Another configuration consists in pipe coiled in circular fashion such that each
loop overlaps the adjacent loop as shown in Figure 2.2-2 (Chiasson, et al. 1999). This
7ı
systems.
Figure 2.2-2: Slinkypipe configuration in a hydronlcally-heated slab
A proper design of a heated slab requires the detennination of the needed heat
flux on the surface based on design weather conditions. To achieve this, it is necessary to
specify the proper spacing ofthe heating elements and the depth at which they are
embedded into the slab. Other parameters such as flow rate, insulation, and operating
temperature also need to be specified to achieve the proper heat flux.
2.3ı Steady State models
In this section, an analysis of the most accepted steady state models is presented.
The steady state snow melting models are either 1-D or 2-D models. The development of
snow melting models started with a classification of service done by Chapman (1956).
Several authors have follow Chapman's work. Some of them based their analysis on I-D
models like Kilkis (1994a, 1994b) and Ramsey et a1. (1999). Other authors developed the
concept further to 2-D models like Schnurr and Rogers (1970). However, even with the
increased use of2-D models (steady state and transient), ASHRAE still presents I-D
8ˇ
steady state model (Ramsey, et al. 1999) for snow melting systems design in the
ASHRAE handbook series.
2.3.1 I-D Steady State Models
These types ofmodels are surface-only models. They are based on the energy
balance at the surface of tbe slab and do not concern themselves with conduction heat
transfer in the slab. The main difference between the different models is the way the heat
losses are evaluated. Chapman (1952, 1956), Kilkis (1994a, 1994b) and Ramsey et al.
(1999) have followed this concept. These models were designed as a guideline for snowmelting
applications with enough accuracy for system design.
Chapman Model
In this model the author originally classifies a heated slab according to the ability
of the system to melt snow. This classification is used to prescribe the operating
conditions of the slab. These classifications are based on the ratio of snow-free area to
total area (Ar) of the slab during the snowfalL The classes are described as follows:
• Class 1. (Residential): The entire surface maybe fully covered with snow
during the snowfall (Ar=O). The system is expected to melt the
accumulated snow some time after the snowfall. Examples of this type are:
residential walks or driveways.
• Class 2. (Commercial): up to 50% of the surface maybe covered with
snow. Examples of this type are: sidewalks, driveways and steps in
hospitals.
9
• Class 3. (Industrial): The entire surface is kept free from snow
accumulation. Examples of this type are: toll plazas in highways,. bridges,
loading areas in airports, etc.
For Ar = 1, the system should melt snow fast enough to prevent any accumulation
over the entire surface. For Ar = 0, the entire surface ofthe slab it is assumed to be
covered with snow thick enough to prevent any heat and evaporative loss. For this reason,
the snow-free area ratio is considered as the main design parameter for snow melting
system designs. This parameter is defiI'\ed as:
(2.1)
Where:
A r = Snow-free area ratio, (dimensionless).
AI = Snow free area, (m2 or :tt2).
As = Snow-covered area, (m2 or fe).
In this model, the correlation used to describe the heat requirement for the snow-melting
process is the steady state energy balance for the required total heat flux at the
upper surface ofa snow-melting slab during snowfall. The correlation used for this
purpose is shown below:
(2.2)
10
2
Where:
qo = total required heat flux per unit area ofthe surface, (BTUIhr.ft2 or W/m ).
qs = Total sensible heat flux, (BTUIhr.ft? or W/m2
).
qm = Heat offusion of snow, (BTUIbr.ft? OJ! W/m2
).
A,. = Snow free area ratio, (non-dimensional).
qh = Sum of convection and radiation losses, (BTUIhr.:tY or W/m2
).
qe = evaporative losses, (BTU/hr.:tY or W/m2
).
The total sensible heat flux required to raise the snow (ice) temperature up to
melting and is calculated as follows:
(2.3)
Where:
Tf= Water film temperature, CF or °C).
Ta = Ambient temperature, COF or °C).
s = rate of snowfall in water equivalent, (in/hr or mm/hr).
The total heat required for a change from solid (snow or ice) to liquid (water)
state can be calculated as:
11
(2.4)
"Where:
C2 = 746 [BTU/(ft2 in)] or 92.65 [W hr /(m2 mm)].
The combined heat flux due to convection and radiation heat transfer over the slab
surface was calculates as:
(2.5)
Where:
Ie = film coefficient ofheat transfer due to conv~tion, [BTU/(h ft2 OF) or
W/(m2 °C)].
fr = film coefficient ofheat transfer for radiation, [BTU/(h WOF) or W/(m2 °C)].
The values for fc and fr are obtained from empirical correlations developed by the
Chapman (1952). These equations were obtained with a wet test slab of4ft? Then the
values obtained under different environmental conditions for wind speed and
temperature- Then. the data was tabulated and plotted. Finally, the least square method
was used to develop the equations. Under these circumstances, the evaporation load could
be overestimated for bigger applications because the small thennal mass of the concrete
slab which was never considered. The equations obtained were:
(2.6)
12
(2.7)
Where:
C3 = 0.27 [BTU/(ft2 mi OF)] or 3.43 [J/(m3 °C)].
C4 = 2.5 [BTU/(hr ft2 OF)] or 14.20 (W/(m2 °C)].
Tp = Surface temperature ofthe panel, (OR or K).
Ta = Ambient temperature, (OR or K).
tp = Surface temperature ofthe panel, (OF or °C).
ta = Ambient temperature, (OF or °C).
v = wind speed velocity, (mph or mls).
The heat flux required to evaporate the water accumulated over the slab surface
can be calculated from an empirical correlation. The derivation of this fonnula followed
the same procedure as in the case fOT convection or radiation heat transfer. After the test
the correlation obtained was:
.(2.8)
Where:
13
Pwv = vapor pressure of water, (inHg or Pa)
Pav = Vapor pressure of moist air, (inHg or Pa)
C6= 0.0201 [BTU/(ft2 mi inHg)] or 4.1 88x 10-5 [J/(m3 Pa)].
C7 = 0.055 [BTU/(hr ft2 inHg)] or 5.134x 10-5 [W/(m2 Pa)].
Chapman (1956) recognized the transient nature of storms. He based all his
calculations on four parameters: velocity ofwind, air temperature, relative humidity of
air, and amount of snowfall. In his study, Chapman (1956) stated, "the true analysis of
heat output depends on an analysis ofthe weather". Moreover, in the determination of
convection and radiation he used a mixed correlation to describe two losses that are
sensitive to two different factors (temperature and wind speed) during different weather
conditions. However, these were the only parameters tested under changing weather
conditions.
The author of this model used collected data from specific cities to develop his
correlations. This is a big limitation even in I-D steady state calculations due to the
impossibility to cover all the regions in a single country and all the possible weather
conditions for each analyzed region. Furthermore, some important parameters such as sky
temperature, sky cover, and relative humidity have not been considered. Finally, there is
no discussion about the heating elements distribution.
14
....._---------'-;------"-'----------------~_._-------- Ki/kis Model
Kilkis (1994a) notes, the minimwn temperature of the system is located halfway
between the heating elements. If the surface temperature is below O°C(32°F) the snow
win not melt or eventually become ice. Therefore, a good design will focus on raising the
effective surface temperature in order to maintain the minimum temperature above O°C
(32°F). As described below, Kilkis uses an analogy to a fin in order to estimate the
minimwn temperature.
The energy balance equation for'the surface of the slab is exactly the same as used
by Chapman (1956) in his model (equation 2-2.). However, the correlations used to
describe radiation, convection and evaporative heat losses are different. KiJkis (1994a)
uses an updated correlation to describe radiation heat loss based on the net long wave
radiation under cloudy skies (equation 2-9) and clear skies (equation 2-10). The
mentioned correlations are:
(2.9)
(2.10)
Where:
C1= 10.3 [BTU/(hr ft?)] or 32.49 (W/m2
).
C2 = 8.14xlO-1O [BTU/(br tt2 °R4
)] or 2.695xlO-8 [W/(m2 K4)].
C3 = 7.316 [BTU/(hr ft2 °R4
}] or 2.423xIO-9 [W/(m2 K4)].
15
C4 - 30.15 [BTU/(hr :ft?)] or 95.11 (W/m2
).
Cs :::; 0.74 [BTU/(hr ~ OR)] or 4.20 [W/(m2 K)J.
Tf= Water fihn temperature, (OR or K)
Ta = Ambient temperature, (OR or K)
In the case of the convection heat losses, the correlation used is based on
experimental data. However~ in his study Kilkis (1994a) used a more conservative
approach and the developed empirical correlation was obtained using a bigger surface
area for the testing slab (16ft?). The correlation used is as follows:.
(2.11)
Where:
C6 = 0.14 [BTU/(mi ft2 OF)] or 1.78 [W/(m3 DC)].
C7 = 0.39 [BTU/(fe OF)] or 2.22 [W/(m2 DC)].
U = conccted wind velocity, (mph or m/s).
Tf= Water fibn temperature, (OP Oli DC).
Ta = Ambient temperature, (OF or DC).
Kilkis (1994a) used a corrected wind speed correlation based on previous
assumptions that the wind speed data is recorded at a different altitude than the altitude 16
the actual application. Therefore, his correlation includes the altitude adjustment factors
as described in AsHRAE Handbook (1989):
(2.12)
Where:
a = Wind speed adjustment factor with respect to height from the ground,
(dimensionless).
Ao = Wind speed adjustment factor with respect to terrain, (dimensionless).
H = Height ofsnow melting surface from ground level, (ft or m).
Href= Wind speed recording height from ground level~ (ft or m).
Urtler = Wind speed from local :meteorological da.t.a, (mph or mls).
The heat flux. required to evaporate the water is caleulatui based 011 required
convective heat flux and the Bowen's ratio:
(213)
Where:
Pal' - Vapor pressure ofmoist air, (inHg or Pa),
P" =Vapor pressure ofsaturated. 'air at film temperature, (inHg ,or Pal.
11
,..........._--~
R = Bowen's ratio, (inHg/°F or Pa/°C).
Tf= Water film temperature, (OF or DC).
Ta = Ambient temperature, COF or DC).
The Bowen's ratio is the relation between lat.ent and sensible heat flux and can be
calculated as: The values for Bowen's ratio (R) and Pal" can be obtained from:
R=CPa
(2.14)
8
Where:
Cg = 9.669xlO-3 (l/OF) or 1.746xl0-2 (l/OC).
Pa = atmospheric pressure, (inHg or Pal.
The values for the vapor pressure of moist in the rage between -28.8 °C (-20°F)
to 1.6 °C (35 OF) air were obtained from the following correlation:
(2.15)
Where:
B= Relative humidity, (dimensionless).
C9= 0.0371 (inHg) or 125.64 CPa).
CIO = I.64xl0-3 (inHgl°F) or 1.197xlOI CPa/DC).
ell = 5.235xlO-S (inHg/°F) or 3.202xl0-1 (Pa/°C).
18
Ta = Ambient temperature, (OF or °C).
The work presented by Chapman (1956) includes weather data from only 10
cities. But, one ofthe main design parameters is the rate of snowfall on the design
location. For this reason, Kilkis (1 994a) developed a statistical correlation to provide a
design tool fOT the rate of snowfall based on available weather data for any location:
s' =[SF cJ .os (2.16)
24 .ow
Where:
SF = Maximum amount of snowfall recorded at a given location in 24 hours, (in
ormm).
C = Class, (dimensionless).
Os = density ofsnow, (lb/ft3 or kg/m\
Ow = density of water, (lb/ft3 or kg/m3
).
In this study the author recognizes that the peak snowfall rate can be as much as
tmee times bigger than a recorded average. Therefore, the data used is the peak snowfall
rate. According to Kilkis's (1994a) study, the peak snowfall rate is the most common data
in many countries. However, the density ofsnow is not C01ll.Inon data in most
meteorological surveys. For this reason, a correlation developed by Adlam (1950) is used
to calculate this parameter:
19
(2.17)
'Where:
Cn = 2.6 (lb/ft?) or 41.65 (kg/m\
CI4 = 0.06 [lb/(ft? OF)] or 1.73[kg/(m3 DC)].
C 2 I5 = 0.0027 [lb/(W °F2)] or 0.1401[kg/(m3 OC )].
Os = density of snow, (lb/ft3 or kg/m3
).
Ta = Ambient temperature, (OF or DC).
The analysis ofthe geometry ofthe slab and the piping distribution was a concern
for Kilkis (1994b). The design temperature is based on the minimum temperature.
However. the calculation of this temperature required a deeper analysis ofthe system
geometry. His approach was the use of a composite fin surrounding the heating element.
In this way Kilkis (l994b), was able to establish the minimum and maximum
temperatures based on geometric infonnation for the slab. A detailed analysis can be
found in Kilkis (1994b) paper. The correlation used is:
(2.18)
T = (Trmx -TJ +7:
min (2.19)
Cos(mW)
Where:
20
n= Composite fin efficiency, (dimensionless).
M = Heating element spacing on centers, (ft or m).
m = Fin coefficient, (ft-I or m-I
).
Ta = Ambient temperature, (OF or °C).
W= Halfof the net spacing between adjacent heating elements, (ft or m).
In Kilkis (1 994b), there is no consideration of side or bottom heat losses. To
account for this type ofenergy loss, an approximation of40% is used to reduce the heat
transferred to the upper surface. Since this approximation is just an estimate on cases
with ground contact and 1 ft dept, a new approach was performed for exposed lower
surface and a new experimental correlation was obtained for this case. The heat lost in the
back of the slab and the required heat flux. to the surface was added to obtain the heat
input required on the source for a certain application. The equation for back loss on
exposed lower surface is:
(2.20)
Where:
Cl6 = 1.14 [BTU/(hr tt2 OF)] or 7.496 [W/(m2oC)].
el7 = 0.13 [BTU/(tt2 mi OF)] or 1.6512 [J/(m3oC)].
Ta = Ambient temperature, (OF or °C).
21
Ts = Surface temperature on the back of the slab if exposed, (OF or °C).
U = corrected wind velocity, (mph or mls).
The experimental validation of this model was never done. The model was
validated against finite-element so utioDS. The results showed an error of 10% compared
with the numerical solution under steady state condition.
Ramsey et at. Model
Ramsey et al. (1999) presented this model; the concept used is that a good snow
melting procedure should be developed as a simple design tool. His work is based on
previous models presented by Chapman (1952) and Killcis (1994a, 1994b). Ramsey used
the same concept in developing the correlations in this model as Kilkis (1994a). "This
roo! must show enough ac:curacy fur engineering adcu!m:ions~ (R.amsey~ et at 1999). As
in the Kilkis (1994:a) model, the onlyd.ifferenoo is the wayin which the heat losses :are
calculated.
In RamseY7 et al (1999) model, the heat losses axe assumed to happen only where
no snow aeoomulation exists over the slab. Any snow-build up over the surface oftbe
slab ads· like insulation for the surface.
The suggested calculation procedure is based on the desiredsnow-ftee area ratio
fur a certain condition as input. By doing this, the heat losses are treated accordingly with
the value selected (Le. ifAr =O.6 is selected, the heat loss will be calculated only over this
area and the heating system adjusted accordingly to this condition). The heating system
22
does not control the condition of the slab. This characteristics and the fact that this model
is a steady state model makes it not useful for simulation purposes.
Heat losses due to convection and radiation are calculated using standard
correlations. Then the heat losses are evaluated together as shown in equation (2.22). In
the convective case, the correlation used corresponds to turbulent convective heat transfer
from an exposed surface under a certain wind speed. This is not necessarily the case in
every operating condition. However, it is a good approximation for exposed surfaces. The
convection heat transfer is calculated as follows:
he =0.037 kZ. JRe~·8 Prl! (2.21)
(
Where:
kair = Thermal conductivity of air at air temperature, [BTU/(h OF) or Wf(rn °C)].
L = Characteristic length ofthe surface, (ft or m)
Rer. = Reynolds number based on characteristic length, (dimensionless).
Pr = Prandtl number of air, taken as Pr=O.7, (dimensionless).
According to Ramsey, et al. (1999), for radiation losses, "the mean radiant
temperature calculations use the same general princip.les as those used in determining the
mean radiant temperature for indoor comfort and heat losses calculations". The
correlation used is shown on equation (2.23):
23
(2.22)
(2.23)
Where:
qh = Sum ofconvection and radiation losses, (BTU/fl? or W/m2
).
he = Convective heat transfer coefficient, [BTIJ/(ft2OR) or W/(m2 K)].
Tf = Liquid film temperature, (R or K).
Ta = Ambient temperature, (R or K).
(1= Stephan-Boltzmann constant = 5.6705xlO-8 [W/(m2K4
)] or 1.7123xlO-9
[BTU/(ft2 R4
)]
es = Emissivity of surface (snow or dry), (dimensionless).
TMR = Mean radiant temperature, (OR or K).
Tcloud = Temperature of the portion of sky that is covered, (OR or K).
Tskyclear = Temperature ofportion of sky that is clear, (OR or K).
Fse = Fraction ofradiation exchange between surface and cloUds, (dimensionless).
A curve fit is used to calculate the temperature of a clear sky. The data
corresponding to the equation fit can be found in Ramsey, et al. (1982). In case ofthe
temperature for the clouds the guidelines of U.S. Standard Atmospheres (1962) was used
24
and the clouds are asswned to be at 3,048 meters (10,000 :ft.). According to these
guidelines, the temperature decreases about 6.4 K per 1000tn (3.5 oR / lOOOft.)..
Therefore, the cloud temperature can be found by subtracting the temperature decrease
from the ambient temperature.
The product ofthe evaporation rate and the heat ofvaporization is used for the
calculation ofthe evaporative heat loss. The evaporation rate can be calculated by
multiplying the mass transfer coefficient with the differential humidity ratio between a
liquid surface film and air. The equations used for this purpose are:
(2.24)
(2.25)
Where:
qe = Evaporative losses, (BTU/m2 or W/m2
)
pa= Density of dry air, (lbJft3 or kga/m3)
hm = Mass transfer coefficient, (ft/s or mls)
Wf= Humidity ratio of saturated air at the film surface temperature,. (lbw/~ba or
Wa = Humidity ratio of ambient air, (lbw/lba or kgwlk&a,)
hfg = Heat of vaporization, (BTUllbwor Jlkgw)
25
Pr = Prantl Number (non-dimensional).
Sc = Schmidt number (non-dimensional)
he = Convective heat transfer coefficient, [BTU/(h fl? OR) or W/(rrC K)]
Cpa = Specific heat capacity ofair, [BTU/(lb OR} or J/(kg K)]
All the parameters used. by these correlations have to be measured using
coincident values for the climate factors (i.e. wind speed, snowfal1rate~ambient
temperature, etc.). This model "only analyzes the upper surface ofthe snow melting
surfac.e" (Ramsey, et aI. 1999). No considerations are taken for edge and back losses
during the analysis. The only criterion shown is that the losses can be considered. in the
range between 4% and 50% depending on the construction design. Moreover, a wide
range ofparameters should be considered. and no correlation with the potential heat loss
is shown.
In the Ramsey, et al (1999) model, no experimental validation was perfonned.
The only validation done is comparison with previous case analyses in six cities. The
error with previous cases ranges from 5% to 15%. Although more cities where included.
on the meteorological survey (38 in total), there is still a great limitation regarding the
information required to perfonn thi.s calculations.
2.3.2 2-D steady state models
On one hand, there are several steady state models that Use a one-dimensional
approach. On the other hand, there are only a few models that take two-dimensional
analysis into consideration (e.g. Schnurr and Rogers 1970, finite-element techniques,
26
etc.). These models require some computational effort to perfonn the calculations.
Nevertheless, these techniques can still be used in snow melting design.
SChnurr and Rogers Model
Schnurr and Rogers's (1970) model was developed using steady state finite
difference equations to describe temperature gradients thru the heated slab. Energy
balance equations were used to calculate the temperature at each nodal point.
Unfortunately, this model uses variable sized grids equal to ~ of the pipe diameter. This
fairly coarse grid needs to be matched only for cases where the distance between pipes is
a multiple ofthe diameter ofthe pipe. Another related problem to the use of a coarse grid
is the impossibility ofmatching the round geometry of the pipes to the square nodal
distribution.
This model does not clearly layout the asswnptions required to model the effects
of the exposed surface. The convection and radiation heat transfer coefficients are
combined in the same way as in the one dimensional steady state model presented by
Chapman (1956). The radiant temperature that is considered in the combined coefficient
is fairly simple and does not take the radiation effects of the surroundings into
consideration. The only considerations made are related to the air temperature and wind
speed.
The heat balance is made following the original model ofChapman (1956). In the
Chapman model simplified correlations were used to evaluate the different heat flux
components. The heat balance equation is shown below:
27
(2.26)
Where:
qo = minimum required heat flux, [BTU/(h.ft?) or W/m2
].
qs = heat flux required for sensible heating of snow, [BTU/(h.ft?) or W/m2
].
qm= heat flux required to melt the ~now, [BTU/(h.ft?) or W/m2
].
qe = heat flux required to vaporize the water film, [Btu/(h.fl?) or W/m2
].
he = combined coefficient for convection and radiation, [Btu/(h.ft?) or W/m2
].
ta = dry bulb air temperature, (OF or °C).
One problem in this model is that it does not consider the snow accumulation over
. the slab surface. Therefore, this cannot be considered as a snow-melting algorithm. All
the analyzed conditions are on surface free of snow during steady state in a c1ass-3 heated
slab. Another problem is that this model only can handle a single homogeneous layer of
material (i.e. concrete). There is no consideration of different properties for materials
used in the slab. The slab is analyzed from the outer surface of the pipe. Finally, the
optimization is valid only in cases where adiabatic conditions are present in the bottom
and sides of the slab. Even so, this model could be adapted to receive a different
boundary condition on the bottom of the slab.
28
Sclmurr and Rogers's (1970) model was able to provide some optimization
compared with previous models. Analysis using this model showed that the surface ofthe
slab is not homogeneous. There is a maximum temperature over the heating element
position and a minimum. temperature in between the pipes, as expected. The model
further showed a strong relationship between the geometry ofthe slab and the surface
temperature distribution. Kilkis (1994b) used this characteristic while developing his
steady state model.
2.4 Transient models
In a steady state model, the dynamic response ofthe slab is not considered. With a
time constant in the order ofhours and since the system works intermittently depending
on the weather conditions, the assumption of a steady state model is too conservative.
Furthennore, the heat to the surface cannot be delivered instantaneously and the transient
effects such as the period between the starting ofthe system and the assumed steady state
can be significant (Rees, et aI., 2002).
Th.e models presented by Chiasson, et al. (1999) and Rees, et al. (2002) give
detailed explanations for how a two-dimensional transient model can accommodate the
pipe geometry and the slab geometry. The first model used a finite difference rectangular
grid. Liu, et at (2002) presented a refined version ofthe Chiasson et al. (1999) model for
the 82nd TRB annual meeting. This second model is based on a finite volume method
with structured boundary fitted grids in order to handle the complex geometries involved.
In any case, all transient models require extensive computational effort. The big
advantage is that transient models can be used for simulation purposes and design
optimization.
29
According to Rees, et al. (2002), a comparison ofthe calculated steady state heat
flux required to maintain the pavement snow-free for a number ofhours and the transient
analysis for the same condition can show that the power required is up to five times
greater than what is projected using a steady state calculation. Another observation shows
that some steady state calculations require ofmore heat flux thru the surface than is
possible to provide in practice. Moreover" accounting for the transient nature ofa storm
can show that storms that start with low loads and then increase in intensity will require
more heat flux than those that start with high loads which and decrease with time.
2.4.1 I-D transient models
Transient swface models have been developed as pan ofthe solution ofmore
detailed models. The only one-dimensional model developed is the one used by Rees, et
al. (2002) to solve during the snow accumulation in their model. This transient model will
be detailed as part ofthe complete 2-0 analysis of their model.
2.4.2 2-D transieBt models
Only few models have managed to account in part for transient conditions. With
the advance of computer science, more models representing a transient behavior have
been developed. Some ofthem are following previous steady state models like Schnurr
and Falk. (1973) that is a natural extension from the previous model (Schnurr and Rogers
(1970)). Unfortunately~in Schnurr and Falk's (1973) model the same assumptions and
boundary conditions from the steady state counterpart are still on use. Since this model
does not allow snow accumulation, it is not considered as a snow-melting model.
However, some models followed the same concept of finite difference analysis (i.e.
Chiasson et aI. (l999) and Uu, et al. (2002», and they do account for snow accumulation
30
and different boundary conditions on the slab. Moreover, new solution techniques are
being developed. Rees, et al. (2002) developed new algorithm. using boundary fitted grids
that allow a detailed analysis ofthe transient behavior ofthe slab and complex geometries
can be bond together.
Leal and Miller Model
Leal and Miller's (1972) model is an extension ofthe Schnurr and Rogers (1970)
steady state model and the first two-dimensional transient modeL The authors ofthis
model developed a transient model based on the same boundary conditions as in the
Chapman model. To accoWlt for the 2-D grid this model uses a polar coordinate system
to generate the equation set. A point matching technique is used to achieve the solution
for these equations. Unfortunately, the authors do not show any inform.ation about the
grid geometry or the solution technique.. For example, inform.ation concerning to the grid
used to accommodate the round tube geometry in a transversally shaped square slab is not
presented.
One ofthe main problems with this model is that the solution is oversimplified.
For example, this model does not consider snow accumulation. Therefore, the model is
not considered a snow-melting model. Moreover, the heated slab should be wet according
to the heat loss equations used (equations (2.3), (2.4), (2.5), (2.8». These equations where
developed experimentally by Chapman (1956) over a wetted. slab and used under this
condition. But, the solution presented by Leal and. Miller (1972) implies that there is no
snow precipitation over the slab and no previous assumptions are made related to a wet
31
surface. If this is the case~ the slab should be dry and not wet, thus making the heat
balance used in the model inappropriate.
Finite Difference Rectangular Grid Model (FD-RG)
Originally, Chiasson et at (1999) presented the FD-RG model. However,. it was
lately refined by Liu, et al. (2002). This model Uses a finite difference in a rectangular
grid with a node-centered approach. This model was implemented in HVACSIM+
software as a component model. The modelnnplementation is somewhat convenient to
use because ofan effective aUI for the end user.
The domain in this model is a section equivalent to one half ofthe pipe spacing.
This is quite reasonable because ofthe symmetry and small temperature differences
between the :flow on adjacent pipes. The y direction corresponds to the thickness ofthe
slab. The x direction corresponds to the distance between the centerline ofa pipe to half
the distance ofthe adjacent pipe. The average top surface temperature obtained in this
domain is considered to be the average surface temperature for the entire slab. This
approximation is made because the model neglects the edge loss of the slab. Figure 2.4-3
shows a representation ofthe model domain. All heat flow is assumed to be into the node.
The governing equation is the two dimensional fonn ofthe transient heat
diffusion equation:
f) 2T a2T 1 aT --+--=-8x2
8);2 a at (2.27)
32
... -- I.. Iı ..j- -·1
lll'Tl8I - HelIt Flux Bol.ndary ConcbJn
Pipe WlIIII
(Fluid Convectlon from irt
pipe flow)
- --,
L.. .I
i
a i
I I
Heat Flux Bou'lcIary - Slab Bottom Surface (convection) r--..
lix.. 1f2 01 pipe spacilg
Figure 2.4-3: Model domain showing finite difference grid lind bOllndary conditions
When this equation is discretized in time, the finite difference form of the Fourier
number appears in all nodal equations. The following equation shows the Fourier number
with square nodal spacing, illl equal to !1y.
alit
Fa =T:":2l (2.28)
\~2 )
Where:
Fa =Fourier number, (dimensionless).
a = Thennal diffusivity, (m2/s or fi2/s).
Llt = Time step, (s).
L1x = Grid size in x direction, (m or ft).
33
This model is not unconditionally stable. As a transient problem, the solutions
must reach the steady-state condition with increasing time. However, according to
Incropera and DeWitt (1996), the transient solution for the explicit finite difference
approach is characterized by induced oscillations. These oscillations are physically
impossible and may become mathematically unstable, causing the solution to diverge
from the actual steady-state conditions. The stability criterion for a 2-D model is given
by:
1
Fa:5:- (2.29)
4
Where:
Fa = Fourier Number, (dimensionless)
If LU is equal to the radius times : ,. the grid domain in the x direction will be a
multiple of the pipe radius and equivalent in length to quarter the perimeter ofthe pipe. In.
overall the sum of the yand x dimensions will be equivalent to the pipe perimeter. To
calculate ~t, the prescribed value of a is evaluated in order to achieve mathematical
stability.
The heat fluxes at the boundary nodes are caused by convection, radiation,
evaporation of rain and melted snow, and melting ofsnow. On the upper surface the
mentioned fluxes are considered according to the temperature and condition of the
surface node. On the lower surface, it is possible to choose between expos€Xl and
perfectly insulated (adiabatic) surfaces. Only convection and radiation heat transfer is
34
calculated at the lower surface if it is exposed. The syuunetry lines found at both sides of
the model domain are by definition zero-flux boundary conditions. The heat flux at the
pipe surface nodes and at the exposed slab surfaces are heat flux. conditions. Inside the
pipe, convection heat transfer causes the heat flux from the fluid. In the upper surface of
the slab, convection, solar radiation, thenna] radiation, sensible and latent heat conditions
must all be accounted for. In the bottom surface ofthe slab only convection is used,
aLthough, if necessary, adiabatic conditions can be set.
The correlations used to calculate the convective heat flux over the surfaces are
shown below:
(2.30)
Where:
q"e = Convective heat flux, [BTU/(h.:ft?) or W/m2
].
he = Convection coefficient, [BTU/(hr fi? OF) or W/(m2 °C)].
Tamb = Ambient temperature, (OF or °C).
T(m,n) = Temperature ofnode (m,n), (OF or °C).
To calculate the convection coefficient a Nusselt number is found. The
appropriate number is calculated using the correlations for natural or forced convection
according to the particular case. For the natural convection case, the correlations used are
those for a heated upper sunace. These correlations are shown below:
35
Nu = 0.54Ra4" I
(104<Ra<107 -laminar flow) (2.31)
Nu = O.15Ral"
1
(1 07<Ra<1011
- turbulent flow) (2.32)
Where:
Nu = Nusselt Number, (dimensionless).
Ra = Rayleigh Number, (dimensionless).
For forced convection heat transfer, ~e correlations used are:
1 I
Nu = 0.664 ReI Pr3 (Laminar flow regime) (2.33)
• I
Nu = 0.037 Re5 PrJ (Mix,ed and turbulent flow regimes) (2.34)
Where:
Nu = Nusselt Number, (dimensionless).
Re = Reynolds Number, (dimensionless).
Pr = Prandtl Number, (dimensionless).
The convection coefficient is then computed by:
Nu·k
hc =-- (2.35)
L
Where:
36
he = Convection coefficient, [BTD/(hr fi? OF) or W/(m2 °C)].
Nu = Nusselt Number, (dimensionless).
k = Thermal conductivity of air at pavement node, [BTU/(hr OF) or W/(m °C)].
L = Characteristic length, (ft or rn).
The radiation heat flux is calculated over the upper and bottom surfaces. The long
wave radiation is calculated using a linearized radiation coefficient shown in equation
(2.36). Then the heat flux is calculated usin&.equation (2.37). If the exposed surface is
the upper one, the temperature used is the sky temperature. If the lower surface is
exposed, an approximation to the ground temperature is made using the air temperature.
(2.36)
(2.37)
Where:
hr = Linearized radiation coefficient, [BTU/(hr ff OR) or W/(m2K)].
c= Emissivity coefficient, (dimensionless).
T(m,n) = Temperature of surface node, COR or K).
T2 = Temperature of sky or air (exposed lower surface), COR. or K). •
37
2
q "r = Radiation heat flux from pavement surface, [BTU/(hr ft2) or W/ro ].
Solar radiation heat flux is considered directly for every time step on the FD-RG
model. This information is directly applied as part ofthe environmental information. This
infonnation is used only over the upper surface ofthe slab.
Evaporative heat fluxes are considered only when the surface is wet or covered
with a mixture ofice and water. This model presumes that there is no rain accumulation
on the surface. Thus, all the rainfall is considered to leave the slab surface immediately
after impact, while still fonning a thin film over the slab surface. In the case ofmelted
snow, the model uses the same approach as specified by Rees, et al. (2002) to
approximate the effect of water being retained in the snow due' to capillarity action.
The evaporative mass flux is calculated by equation (2.38) at each surface node:
(2.38)
Where:
Wair = Humidity ratio of the ambient air, (dimensionless).
W(m.l) = Humidity ratio of the saturated air at a top surface node, (dimensionless).
hd = Mass transfer coefficient, [lb/(ft2 s) or kg/(m2s)].
The mass transfer coefficient is calcu~ated with the Chilton-Colburn analogy by
the following equation:
38
(2.39)
Where:
he = Convection heat transfer coefficient, [BTU/(h ft? OF) or W/(m2 °C)].
Cp = Specific heat capacity of air at node temperature, [BTU/(lb OF) or
J/(kg °C)].
Le = Lewis number, (dimensionless)..
The Lewis number is the ratio between thennal diffusivity (n) and the Binary
mass diffusion coefficient (Dab). This last value is calculated using the correlations
Marrero and Masson (1972):
CT2.072
D =---!...I_- (2.40)
ab P.
OJT
Where:
T = Absolute temperature ofair, (OR or K).
Pair = Pressure of air, (inHg or Pa).
The evaporative heat flux due to evaporation is calculated as follows:
qe"v ap = hfg 'm"evap (2.41)
39
Where:
hfg = heat of evaporation, (BTUlib or Jlkg).
riz':vap = Evaporative mass flux, [lb/(ft? hr) or kg/(m2 s»).
The last heat flux value evaluated at the top surface corresponds to the
combination of sensible and fusion heat required to melt the snow. The sensible heat flux
corresponds to the energy required to elevate the temperature of tbe snow to ooe (32°).
This value is calculated as follows:
(2.42)
Where:
rh;re = Mass flux of water equivalent ofthe snow or freezing rain, [lb/(ft2 s) or
kg/(m2 s)].
C pw = Specific heat of water, [BTU/(lb OF) or J/(kg °C)]
T air = Temperature of ambient air, (OF or °C)
T(m, I) = Temperature ofsurface node, (OF or °C).
The fusion heat flux is the flux required to produce the change ofstate for the ice.
It is determined with a heat and mass balance on a specified top surface node. The model
evaluates the actual mass flow rate at each time step. This value is calculated as the
40
minimum between the potential snowmelt flow rate and the accumulated mass of ice on
each time step. This equation is shown below:
." " 1/ mice-accumulated-current (2.43) mmelt-acluil/ = min(mmell-patellrlal' " 111 J
." (." 0) (2.44) mmeJt.-pateTltial =max mmell_heaIBa/ance.
n " " """ • " qsolar + qTad + qcony + qsell + qevap + qcond.SlFif (2.45) mmell healBalance =
Where:
q;olar = Heat flux due to solar radiation, [BTU/(h ft?) or W/rn2
].
q;ad = Heat flux due to long-wave radiation, [BTU/(h ft2) or W/m2].
q:OllV = Heat flux due to convection, [BTU/(h fl?) or W/m2
].
q;en = Sensible heat flux on snow, [BTU/(h ft2) or W/m2].
q';.,op = Evaporative heat flux of slush or water, [BTU/(h ft2) or W/m2].
q:nd.,suif = Heat flux conducted to the surface of the slab, [BTU/(h ft2) or W/m2].
hif= Latent heat of fusion of water, (BTU/lb or J/kg).
This model is designed for simulation only ofhydronic systems. However, the
model could be modified to account directly for the heat flux ofa buried electric cable.
41
At present time, this boundary heat flux cannot be passed directly to the slab. This
boundary condition is applied in the nodes where the internal piping is present.
In order to model heat flux due to heat exchange from the fluid the following
correlation is applied:
(2.46)
Where:
q "fluid = Heat flux through the pipe wail, [BTU/(h ft2) or W/m2].
U = Overall heat transfer coefficient, [BTIJ/(h ttl OF) or W/(m2 °C)].
T(x.yj = Non-surface node temperature, (OF or °C).
The overall heat transfer coeffici~':mt can be expressed as:
1
U=---- (2.47)
_1_+_1_
hp ..id kp1p<
Where:
U = Overall heat transfer coefficient, [W/(m2 °C)].
h fluid = Convection heat transfer coefficient for fluid, [BTU/(h ft? OF) or
W/(m2oC)].
k pipe= Thermal conductivity ofpipe material, [BTU/(h ft OF) or W/(m °C)].
I = Length, (ft or m)..
42
The convection heat transfer coefficient is given by the following relationship:
Nu·k h = fl (2.48) c
L
Where:
he = Convection heat transfer coefficient for fluid, [BTU/(h ft? OF) or W/(m
2
°C)].
kiF Thermal conductivity of fluid, [BTU/(h fe OF) or W/{m2 °C)].
L = characteristic length (defined as the inner diameter of the pipe), (ft or m).
According to Liu, et a1. (2002), the Nusselt Number (Nu) is com.puted with the
Gnidinki equation as shown below:
(2.49)
Where:
Nu = Nusselt Number, (dimensionless).
f= Friction factor, (dimensionless).
Re = Reynolds Number, (dimensionless).
Pr = Prandtl Number, (dimensionless).
The fiicti.on factor is given by:
43
f = [1.58ln(Re)-3.28]-2 (2.50)
Where:
Re = Reynolds Number, (dimensionless).
The fluids that can be used are water or antifreeze solutions. The thennal
properties for the working fluid are calculated for each time step. The subroutine that
performs this computation is based on the ASHRAE Handbook of fundamentals (SI)
1997. The !hennal properties for water, Propylene Glycol and Ethylene Glycol can be
calculated for a given average fluid temperature and antifreeze (% volume). The average
fluid temperature is calculated iteratively using as initial value the previous time step
converged value.
Finite Volume - Boundary Fitted Grid Model (FV-BFG)
In 1999 the ASHRAE project 926-RP entitled "Snow Melting Algorithms and
Data for Locations Around the World" was presented. The objective of this project was
to update the available data for the HVAC Applications Handbook (ASHRAE 1995). The
ASHRAE 1090-RP is an extension of the 926-RP. The ASHRAE 926-RP presents a
steady state load calculation procedure (Ramsey, et al. (1999) model) in order to find an
instantaneous flux required to produce a given snow free area ratio. ASHRAE 1090-RP
was developed by Rees, et al. (2002) and uses as many of the correlations previously
developed as possible. However, one of the main objectives ofthis current project is to
develop a two-dimensional transient analysis simulation tool. Under these conditions, the
44
free area ratio should be treated as an output for the model and not as an input,. as is in the
ASHRAE 926-RP.
Three parts fonn this model; the user interface, the finite volume solver, and the
boundary condition mOdel. The user interface is used to pass all the required information
to the finite volume solver and the boundary condition model. The finite volume solver
has two components, the grid generation module and the surface temperature calculator.
The boundary condition Inodel is used to evaluate the condition ofthe slab and provide
information to the finite volume solver to calc,u late the surface temperature. This last
component is considered as the I-D transient model associated to the finite volume solver
used to find iteratively a solution for the surface temperature. .
The finite volume solver model is based on the integral fonn ofthe partial
differential equation ofthe Fourier equation for heat conduction, as shown below:
~ f¢dV = Jrv¢ndS (2.51 )
at v s
Where:
¢= Temperature.
r = Thermal diffusivity.
v = Control volU1Ile.
S = Surface ofthe control volume.
n = vector nonna! tothe surface ofthe control volume.
45
The right hand side of the equation represents the diffUsion fluxes. In its discrete
form the sum ofthe diffusion fluxes through each cell face can be repres,ented as:
Jrv¢ndS:::: LF;D (2.52)
S i
Where:
v ¢ = Temperature gradient
S = Surface of the control volume
n = vector Donnal to the surface of the control vohnne
i = north, south, west, east (for a 2-D cell)
D = Count for the diffusion term
On a particular cell the equation (2.52) is transformed into:
(2.53)
Where:
Se = Area of the east face
The main challenge in using equation (2.53) is the calculation of the temperature
gradient (Vrji) at each cell face. To overcome this difficulty, Rees, et al. (2002) use the
values ofthe variable at the cell centroid ((Jp and fA) and the distance between these
46
points (Lp'E'). This, is a valid approximation if the cell is orthogonal. The equation (2.53)
is reduced to:
F D =r s (8¢J (2.54) e ee '
8q e'
As expressed by Rees, et aI. (2002), the main idea is to "preserve the second order
accuracy by making the calculation of the gradient along the nonnal face". The use ofthe
values ofthe variable at the p' and E' is essential. See Figure 2.4-4.
,ne
E'
""' e' ~ IE
se r
YL
x
Figure 2.4-4: The local coordinate system at the eastface ofatypicalflnitevolume cell
The values at these points are calculated using a "deferred correction" approach,
as follows:
(2.55)
47
This is done using the previous values ofthe variable ¢. When the solutions approach
convergence, the gradient tenns along qcancel and the gradient along the normal to the
face is preserved as desired. The (o«jY&)n term can be explicitly calculated from the
central difference as:
(2.56)
Where:
LrE ' = Distance between p' and E'
t/Jp.;t/JE. = Variable evaluated at P'; E'
To obtain the finite volume solution, the partial differential equation needs to be
integrated with respect to time. Rees, et a1. (2002) used a first-order backward
differencing approach in an implicit formulation. This procedure is fully described by
Rees (2000). This approach results in an unconditionally stable equation. It is said to be
first-order accurate on time and second-order accurate on space. This equation is shown
below:
(2.57)
Where:
n = Variables at previous time step
48
n+1 ::;:: Variables at current time step
~t;;;;; time step
N,S,W,E:;:: sub index used to denote the face ofthe cell.
The final equation for a control volume can be determined ailer the integration
and discretization procedures have been applied. Each control volume will have the fonn:
.(? 5'S·) ' : ".......n
Where:
a , b = coefficients ofthe algebraic equations
rPnb = Variable at defined cell
n = North, south, west, east
The ap value for one control volume becomes Clnb for the next cell. Therefore, in a
two dimensional model, a penta-diagonal matrix is formed and can be solved with a
convenient matrix solver technique. The method used by Ress, et a1. (2002) is the
Strongly Implicit Method (Stone, 1968).
The grid is generated using blocks of cells. The blocks are described as a regular
rectangular array of cells. Each interconnected block has the same number of cells. Thus,
the shape of the edge (i.e. lines, arcs) is not particularly important. In general, it is
necessary to have at least 4 blocks to represent the pipe and pavement layer CO!l.taining
the pipe. More blocks can be used ifmore layers are necessary. See figure below:
49
n+1 = Variables at current time step
~t = time step
N,S,W,E = sub index used to denote the face of the cell.
The final equation for a control volume can be determined after the integration
and discretization procedures have been applied. Each control volume will have the form:
(2.58)
II
Where:
a, b = coefficients ofthe algebraic equations
tAb = Variable at defined cell
n = North, south, west, east
The ap value for one control volume becomes anb for the next cell. Therefore, in a
two dimensional model, a penta-diagonal matrix is formed and can be solved with a
convenient matrix solver technique. The method used by Ress, et a1. (2002) is the
Strongly Implicit Method (Stone, 1968).
The grid is generated using blocks of cells. The blocks are described as a regular
rectangular array of cells. Each interconnected block has the SaIne nUmber of cells. Thus,
the shape of the edge (i.e. lines, arcs) is not particularly important. In general, it is
necessary to have at least 4 blocks to represent the pipe and pavement layer containing
the pipe. More blocks can be used ifmore layers are necessary. See figure below:
49
2
3
"
5
Figure 2.4-5: Block definitions on a typical grid layout
The cell distribution on the grid must be optimized. To optimize the grid a higher
cell density is applied where higher temperature gradients are expected. Other regions use
large cell sizes. The density of the cells determines the accuracy of the solution.
However, there is also a strong relation between the density of the cells and the
computational time required. The higher the ceU density, the more computational time
required.
According to Spitler, et aL (2001), a higher density grid for this application is
desirable close to the pipe. It is in this pomt where the higher temperature gradients are
found. To follow this criterion, an exponential distribution of the grid is applied between
the corresponding block edges. The equation used for the distribution is shown below:
Le:x:poMru
e-n - -1
exponential. = ---- (2.59)
I e exponent 1 50
Where~
i = cell counter that varies from i to n-l (number ofceUs in distribution line).
exponent = normalized distribution along the line or arc.
The grid generation rests on several prior detenninations.In order to ensure there
is enough space for a correct grid generation the value ofpipe spacing and depth under
the surface are checked. The minimum value can be shown to be 1.5 times the outer
diameter of the pipe. An example of th.e grid geperated for this model is shown in. Figure
2.4-6.
The estimated number of cells is obtained from the geoIIletric description ofthe
slab. The model determines the number of cells required for slahs ofdifferent sizes. A
full description ofthe cell generation procedure can be obtained fonn Spitler, et al.
(2001).
51
...-..,JIII/Il
E;:k.L 1I11 J
tf\ \\\\
H-l \ \\\\ \\
Figure 2.4-6: Grid generated/or hydronic system
The finite volume solver is coupled with the boundary condition models by
exchanging infOlm.ation about the surface temperature. However, since the temperature is
constant at the melting point it is necessary to use the heat flux calculated by the finite
volume solver. The boundary condition solver is used for the surface temperature
calculation. This process is iteratively repeated until the temperature converges. At
convergence, the heat flux calculated by the finite volume solver is consistent with the
temperature calculated by the boundary condition model.
The boundary conditions used in this model are similar to those used in the
previously analyzed model (FD-RG). The boundary condition model is a collection of
models. The heat balance is satisfied at the surface nodes according to the cond'tion of
~
the slab. The conditions on the surface ofthe slab can be: dry, wet, dry snow, slush, snow
52
and slush, solid ice, and solid ice and water. These conditions are summarized in Table
2.4-1. Probably the most complex and interesting model is the model describing snow
melting.
Surface condition Properties
Dry Surface free of liquid.
Wet Surface with liquid and above the freezing point.
Dry Snow Fresh fallen snow with no liquid or previous precipitation.
Surface below the freezing point so no snow melting occurs.
Slush Surface contains ice with snow crystals that are fully saturated
with water. Surface temperature at freezing point.
Surface contains snow that is partially melted. The lower part of
Snow and slush snow is saturated witli water and the upper part is fresh snow.
Surface temperature at freezing point.
Frozen water over the surface of the slab. Surface temperature
Solid ice
under the freezing point.
This condition can occur when the ice is being melted or if there
Solid ice and water
has been a rain even over the icy condition.
Table 2.4-1: Surface conditions over the heated slab
Heat and II'lass transfer is calculated for each cell ofthe grid on each time step.
According to Rees~ et al. (2002) the calculation of the height ofthe saturated layer in the
snow pennits calculation ofthe mass of dry snow. So, in addition to heat balance, it is
necessary to keep track of mass transfer for each time step. In order to keep track of the
surface condition the slab is considered dry unless any of the conditions that would
change this dry condition arise.
The 1-» snow-melting model
As described by Rees, et aL (2002), during the melting process the snow can be
considered as a layer of "dry" snow (only crystals with no liquid) and a layer ofsaturated
53
snow (slush) close to the slab surface. Both are considered to be porous nIe~ one with
air and the other with water in the void space between crystals.
Mass transfer occurs as shown in Figure 2.4-7. The quantity of snoW and rain is
detennined from weather data. Snowmelt rate is detennined by energy balaJlce. As a
result, the slush line moves such that previously dry snow becomes slush. So:rne mass is
exchanged with the snow layer, and some mass change results from evaporation ofthe
liquid layer.
Snowfall
Rairdall
Atlnosphere
Snow l<wer
.....f---- Runoff
Figure ,2.4-7: Mass transfer on snow Melthrgmodel
A three-node model is used to model the heat transfer for the- snow-melting
mOdel. The first node is located at the upper surface ofthe snow layer. The second node
is located at the center ofthe snow layer. The last one is located at the saturat'ed slush
layet. This model is represented, in Figure 2.4~.
54
Atmosphere
Snow layer
-+-,-----------~
..c::; : I Saturated (slush) layer
Slab
Conduction
Figure 2.4-8: Schematic representation ofheat transfer in the three·node snowmelt model
The assumptions of this model are:
• Unifonn temperature in the slusb/liquid layer.
• Melting of snow occurs at the lower node only.
• Transfer of solid snow fonn the snow layer to the slush layer is explicitly
accounted by the mass balance.
• While convection from the upper surface ofthe snow is accounted for,
convection due to airflow through the porous.snow layer is neglected.
(The model can include this condition if necessary).
• Convection and evaporation from the slush layer are neglected when
covered with a layer ofdry snow.
55
• Rainfall occurring after a snow layer has fonned is accounted for directly
only at the saturated layer.
• The snow melting process is treated as a quasi-one dimensional process.
This model is fanned by five primary equations: mass balance for the solid ice, a
mass balance for liquid water, and a heat balance on each node.
Mass balance for ice:
dmice " If ---m -mmel/ (2.60)
- s'lOwfall dB
Where:
mice = Mass of snow per unit area in the snow layer, (kg/m2 or Ib/ft2
)
B= Time, (8 or hr)
m;',owfall = Snowfall rate in mass per unit area, [kg/(m2s) or Ib/(ft2h)]
m:e1t = Rate of the snow transferred to the slush in solid form., [kg/(m2s) or
lb/(ft?h)]
Mass balance on the liquid is given by:
(2.61)
Where:
56
() = Time, (s) or (h)
m;aln = Rainfall rate in mass per unit area, [kg/(m2s) or Ib/(tt2h)]
m:W1 = Snowmelt rate in mass per unit area, [kg/(m2s) or Ib/(tt2h)]
m':n"oJf = Rate ofrunoff in mass per unit area, [kg/(m2s) or Ib/(tt2h)]
The amount of runoff is defined as the excess water on the slush layer due to
melting. To take account for the capillarity action, the runoff is limited to 10% of the melt
rate. This consideration is done until the saturated layer has 5.08 em (2 in.). When this
condition is reached the runoff is increased to the melt rate to prevent more water to be
retained. This consideration is done due to the fact that capillary forces reach balance
with gravitational forces. This assumptions where experimentally verified by
Hockersmith (2002).
The total mass of snow and slush are required for the heat balance. This is done
using a relative density of the snow and calculating the thickness ofthe layers. The total
height is calculated by:
h mice (2.62) 10101 - ( )
Pice 1-1]eJf
Where:
h/olal = The total thickness of the snow and saturated layers, (m or ft).
57
mice = Mass ofsnow per unit area in the snow layer. (kg/m2 or lb/tY)
TJeff= The effective porosity of the ice matrix (applies to both layers),
(dimensionless).
Pice= The density of ice. (kg/m3 or lb/ftJ
)
In the same way the thickness ofthe liquid layer is calculated from the mass of
liquid:
(2.63)
Where:
htotal = The total thickness of saturated layers. (m or ft)
m I = Mass of liquid water per unit area in the slush layer, (kg/m2 or Ib/fY)
PI = density ofsaturated layer, (kg/m3 or Ibm/ft3
)
neff= The effective porosity ofthe ice matrix, (dimensionless)
The height of snow is calculated by subtracting the height of saturated layer from
the total height. Once this value is found the mass of snow can be calculated. Thus, it is
used on the energy balance equations.
m dTsnow - " -" -" " snow CP dB - qconduction ,snow qsnowfall qconvection - qradiation (2.64)
58
The conduction heat flux ofthe slush layer to the snow layer is calculated as:
q"cO llduction.snOMl = ksnow (Tslush - Tsnow ) (2.65)
O.5h snow
The heat flux due to new snow over the surface is given by:
(2.66)
The convective heat flux in the upper surface of the snow is given by:
(2.67)
The radiative heat flux is given by:
(2.68)
Where:
ksnow
= conductivity of snow, [W/(m K) or Btu/(h ft OF)]
hsnow = height of snow, (m or ft)
Tsluh = Temperature of slush, (Oe or OF)
Tsnow = Temperature of snow, (OC or OF)
Ta
= Temperature of air, (Oe or OF)
Tsurftce = Temperature of surfaoe, cae or OF)
59
TMR = Mean radiant temperature, (OC or OF)
es = Emissivity ofsurface (snow or dry), (dimensionless)
he =heat convection coefficient, [W/(m2DC) or Btu/(h:fY OF)]. (Convective heat
transfer coefficient is calculated as described in equations (2.31) to (2.35»)
An iterative process is used to calculate the snow surface temperature. Since the
surface temperature is used in the convection and radiation heat calculations. the
equations (2.67), (2.68), and (2.69) are used for this purpose. Then, the energy balance of
the slush node is used since at this point the mixture ofliquid ice -is at equilibrium. The
energy balance is shown in equation (2.70).
T T O.5hsnow ( " ")
sllrface = snow - qconvectiorr + qrodio./ion (2.69)
ksnow
Energy balance at slush node:
m"melJ hIf ="qco ndllctwn ,slab + q"ra inf all "- q con dllctlon,sllow (2.70)
Where:
m:,elt = Snowmelt rate in mass per unit area, [kglm2s or lb/ft?h]
hif= Latent heat of fusion ofwater, [Jlkg or Btu/(h.1b)]
Assuming that rainfall is at air temperature, the heat flux due to rain is given by:
60
(2.71)
Where:
q "raillfafl= Heat flux due to rain, [W/m2 or BTU/(hr ft2)]
m"rainfall = Rainfall in mass per unit area, [kg/m2s or Ib/ft2h]
Cp water = Specific heat ofwater, [J/(kg °C) or BTU/{lb OF)]
Ta = Temperature ofair, (OC or OF)
Tsluh = Temperature of slush, (OC or OF)
This model has been only partially tested. Hockersmith (2002) used the snowmelting
correlation ofthis model for his experimental investigation. However, his
experimental apparatus was a lab-scale flat electrical heater. Therefore, the numerical
solver used to calculate the heat flux on the surface was not used. In Hockersmith's
(2002) investigation, the power used by the heater was directly passed to the snowmelting
correlations. Thus, only the snow-melting correlations where tested.
2.5 Experimental investigations
There are several experimental investigation perfonned. Mainly, experim,ental
investigations of snow melting have been done in cases where industrial de-icing or snow
melting is required over industrial equipment using direct steam or hot air blowers in
experimental designs. In the case of snow melting applications on slabs limited
experimental validation is perfonned. The authors ofthe snow melting models do most of
61
the experimental investigation while developing their models. However, two
experimental analyses are considered in this study.
The first one is the experimental validation presented by Aoki,. et a!. (1987). The
paper presented showed the investigation ofsnow melting process by heating a plate
covered with snow. The second study is a thesis presented by Hockersmith (2002). In this
report an experimental and computational analysis ofthe I-D transient snow-melting
algorithm used by Rees, et a1. (2002) in their model is discussed.
2.5.1 Snow melting by heating from the bottom (Aok4 etal. 1987)
In this experimental investigation the researchers based ~eir work on the concept
ofcapillarity rise ofa fluid on porous media. This rise is caused by the difference in
pressure ofthe hydrostatic column ofthe liquid with the one produced by the surface
tension ofthe liquid in contact with the porous media inside the void space. This is the
same phenomenon that occurs while snow is melting. The capillarity pressure is balanced
with the hydrostatic pressure of the melted water. When the water colunm reaches the
equilibrium point, it will overspiU outside the snow media. Aoki uses this concept to
formulate his model based on parallel capillary tubes of different sizes packed together.
The researchers heat up a plate with snow on it and observe the different melting
process occurring on the snow. By doing these experiments they classified the snow in 3
•
classes depending on the air temperature surrounding the experimental snow:
62
Class A
Occurs when the temperature ofthe surroundings is near DOC (32°F). The melting
process occurs after the plate has been heated up. While the snow melts, the water rises
due to capillary action of the porous media. The water infiltration moves upward. If the
snowmelt rate is bigger than the water penneation the water drainage occurs. There are
two subtypes in this class. Al occurs when the snowmelt is bigger than the water
permeation. Conversely, A2 occurs when the waterfront reaches the upper surface of the
snow.
ClassB
This type ofsnow melting is produced when the temperature of the surroundings
is less than the freezing point. In this case are-freezing of the melted water occurs. After
the plate is heated, the snow melts and the water penneates by capillary action. During
the time when the water raises a portion of the snow layer freezes. The freezing does
decrease the penneation ofwater. At this point, two subclasses are recognized. Bl occurs
ifthe water penneation front has reached the upper surface of snow. In this case, an ice
lens is formed over the snow surface. B2 occurs if the drainage happens before the water
penn.eation front reaches the upper snow surface
Class C
This class is recognized in the case that the waterfront freezes inside the snow.
The permeation rate is lower than the freezing rate. At this point the waterfront decreases
and until a point where the freezing rate is lower than the penneation rate. Wheu this
63
condition has arisen, the problem can become a class "B" or "e" again until all the snow
is melted.
Aoki, et al. (1987) based its analysis on snowmelt Class "C". The equations llsed
to describe this class are divided according to the different layers of the model. The first
layer is the water permeation layer where the temperature is considered constant at O°C
(32°F). In this layer a transient mass balance is perfonned:
dSwp dm dS
ernw --=---ux 8r.'"w ~ (2.72)
dt dx dx
Where:
&= Porosity, (dimensionless)
Pw = Density of water, (kg/m3 or Ib/ft3
)
x = Height, (m or ft)
2 m = Water permeation flux, [kg/(hr m2
) or lb/(hr m )]
Ux = Velocity due to the reduction in volume, (mIhr or ft/hr)
Swp = Water content ratio, (dimensionless)
The second layer is the ice layer, where the temperamrechanges with time.
Convection and conduction effects are considered.
dI; d2I; dI; -=a;---u,,(
2.73)
dt dx2 dx
64
Where:
7;. = Temperature of ice layer, (Oe or OF)
t = time, (hr)
ai = Thermal diffusivity of ice, (m2/hr or fPIhr)
Ux = velocity due to the reduction in volume, (mIhr or ft/hr)
x = Height, (m or ft)
The third and last layer is the "virgin snow" layer. In this layer, the temperature
changes on time. TIris layer also includes convective and conductive effects but in this
case a porous matrix is used.
2
dTp =a d Tp -u dTp (2.74)
dt p dx2 P dx
Where:
Tp = Temperature ofporous layer, (Oe or OF)
t = time, (hr)
ap = Thennal diffusivity ofthe porous layer, (m2/hr)
Ux = velocity due to the reduction in volume, (rnIhr)
x = Height, (m or ft)
65
upper surface ofthe snow and constant heat flux and 100% water saturation on the
bottom surface in contact with the plate. The equations are solved with a finite difference
method and variable space network to move the boundary layers with the layers.
This study shows that heat loss from the upper surface vary depending on the
snow pattern. Therefore, the melting time required by the snow depends on the water
permeation and the saturated layer refreezing. The model was also able to estimate the
water drained from each layer with the experimental observation.
2.5.2ı Hockersmith experimental investigation
The experimentation performed by Hockersmith (2002), is based on an extensive
research ofthe snow melting theory. The properties ofsnow and its metamorphism are
taking into consideration on this document. Metamorphism is a characteristic of snow
due to the change in shape ofthe snow crystals during the melting process.
For his experiments, Hockersmith (2002) used an electrically heated plate to melt
snow from the bottom in a controlled environmental chamber. Measurements on the
crystal size, water runoff, capillarity height of snow (maximum saturation layer height),
and plate temperature under different values ofheat flux were evaluated. The results from
these experiments were compared with the 1-D transient model used by Rees, et al.
(2002).
The results ofthe experimentation showed that the model predict a linear snow
melting rate. However, compared to the experimental data, the snow-melting rate
decelerated during the last hour ofthe experiment. Hockersmith (2002) attribute this
66ˇ
with the snow metamorphism when the water has saturated the snow (during the slush
period). Nevertheless,. the 15% difference observed during the experimentation was also
associated with uncertainties due to the model itself. In the model, side and bottom heat
loss were not considered.
Another important observation was the water runoff model. This parameter has a
direct influence on the maximum saturation layer height. When compared with the
experimentation, a time delay between the modeled runoff and the experimental runoff
was confirmed with an error ofup to 20%. To correct this problem, Hockersmith (2002),
include a simple runoff model that drained water form the slab at a constant rate. The
error then was reduced to 5% for the same case scenario.
67ı
CHAPTER 3
II!
EXPE~NTALAPPARATUS
i
f
3.1 Introduction
The experimental apparatus consists primarily ofa hydronically heated concrete
slab and a snow-making environmental chamber. The concrete slab was embedded with
hydronictubing and is connected to a circulating pump and a controllable electric heater.
The slab is instrumented with a number ofthermocouples. Power input to the heater and
the pump is measured as function of flow rate and differential temperature ofthe fluid.
The environmental chamber is used to control the conditions surrounding the slab.
In the environmental chamber, it is possible to maintain temperatures around the concrete
slab and to simulate snowfall over its surface. Air and surface temperatures were
measured at a number oflocations in the chamber. A detailed description of each
component is given in this section. A data logger is used to collect the infonnation
provided by the slab and from the environmental chamber.
Two auxiliary experimental apparati were required. The first One was used to
measure the conductivity ofthe concrete used in the slab. The second one was used to
measure the specific heat of the concrete. These two important parameters are basic data
68
required by the snow melting models. The experimental apparatus used in each case are
described at the end of this section.
3.2 Concrete Slab
The concrete slab was constructed using a wooden frame (mold), which
subsequently remains attached to the slab. The plywood is 0.0191 m. (0.75 in) thick. The
void space is filled up with fast setting concrete (Quikrete). The dimensions ofthe
concrete are: 0.914 m (36 in) length, 0.864 m (34 in) wide, and 0.152 rn. (6 in) depth. The
length dimension is said to be the same direction as the hydronic tubing.
Three paraUe1lines ofhydronic tubing were located inside the frame in the
direction of the length as shown in Figure 3.2-1. The distance between the pipes is 29 em.
(11.5 in) taking as reference the pipe located in the center ofthe frame. The pipe used
has an outer diameter of 15.87mm (5/8 in) and wall thickness of 2.4 mm. (3/32 in). The
specification for the pipe material is Polyethylene PEX, SDR-9.
69
Figure 3.2-1: Thennocoupk and tubing placement during concrete'pounng process
Thermocouples are placed in several positions across the central section ofthe
concrete slab as shown in Figure 3.2-1. Thermocouple positions are given below in
section 3.2.3.
Concrete was poured inside the wooden frame to create the concrete slab. The
tubes were held in position during the construction process by rigid aluminum pipes
inserted into the tubing. The thermocouples were hold in position using taut line. The
freshly poured product is shown in Figure 3.2-2.
/
70
Figure 3.2-1: Thermocouple and tubingplacement during concrete pouring process
Thennocouples are placed in several positions across the central section of the
concrete slab as shown in Figure 3.2-1. Thennocouple positions are given below in
section 3.2.3.
Concrete was poured inside the wooden frame to create the concrete slab. The
tubes were held in position during the construction process by rigid aluminum pipes
inserted into the tubing. The thermocouples were hold in position using taut line. The
freshly poured product is shown in Figure 3.2-2.
70
3.2.1 Base
An insulated base on casters supports the concrete slab. This section was
constructed with 0.019 m (0.75 in) thick plywood. This section is filled with 0.254 m (10
in) of extruded polystyrene. Plexiglas was used under the central section where the slab is
placed in order to prevent that concrete from flowing into the insulation during the
construction process. The total dimensions ofthe base are 0.851 m (33.5 in) wide, 1.251
m (49.25 in) long, and 0.279 m (11 in) tall.
The base and slab are covered on the sides .by O.127m (5 in) of extruded
polystyrene, which is covered with 0.019 m (0.75 in) thick plywood. A diagram of the
slab and base is shown in Figure 3.2-4. The side covers are held on,tightly by rubber
grips and sealed with silicon in ord.er to prevent circulation of air between the side cover
and the slab.
r SIDE COVERS
CONCRETE SLAB
I
I • ~ _ _ .L ..L _
LOASE LOASE
SIDE VIEW iFRONTVIEW
Figure 3.2-4: Concrete slab and base
72
The whole system needs to be transportable in order to move it in and out of the
controlled chamber. Therefore, two pairs ofwheels are attached to the base. The first pair
is located in the 0.152 m (6 in) from the front side ofthe base. This pair of wheels is
0.254 m (10 in) in diameter. They are attached to an axle, which is fixed to the bottom of
the base. The second pair is located at 0.1016 m (4 in) from the backside of the base. A
detailed central section for the base construction is shown in appendix C.
3.2.2 Heating System
The hydromc heater included in the slab uses Ethylene Glycol (50% volumetric
concentration) as working fluid. It is circulated using a Grundfos pump, model 2P41 o. An
AC electric heater was also integrated in the piping system and allows adjustment of the
power input to the working fluid. The maximwn capacity of the heater is 1500 W. In
order to measure the flow rate in the system, an in-line liquid flow meter was used.
Additionally, two calibrated thermocouples were used. One was placed in the inlet pipe
ofthe system. The second is positioned on the outlet pipe of the system. A simplified
piping diagram of the system is shown in Figure 3.2-5.
73
Concrete Slab
J Piping
r - - - --
1
I
I I ~ Heater
I J To purge tank
I
I
I
I
L - -
I
I
I
I
I
_..J From Purge tank
CD
@
Flow Meter
Recirculation Pump
-t>J<r Valve
'2) Temperature
Figure 3.2-5: Concrete Slab simplifiedpiping diagram
Before the fIrst run of the system it is necessary to remove tramped air from the
system. For this reason, a set ofball valves was used before the inlet ofthe pump.. The
ball valves allow the user to connect a purge tank where the working fluid can be
separated from the air and circulated. After the first run the purge tank is removed and the
valves aligned for flow through the slab.
The free space in the front of the concrete slab is used to place all the auxil'ary
equipment. The auxiliary equipment includes the circulation pwnp, the in-line liquid flow
tneter, the heater, connection piping, and in-flow thennocouples. The auxiliary equipment
was covered with insulation. However, it was necessary to reduce the insulation for the
slab in this point in order to make room for the auxiliary equipment.
Since the main elements (i.e. heater, pump) require almost all the space available
it is not possible to cover tIlls side with a wooden cover. As a result, the insulation around
these elements is only 38.1 rom (I.5in) thick and additional energy loss is expected from
74
with insulation.
Figure 3.2~: 8) AlIXiliary equipment, b) fUlXiJiary equipment covered with insulation ;nsUh the snow
chamber
The pump and the heater produce the power input in the system. The power input
produced by the pump is regulated with a VARIAC (Variable ratio transformer) device.
The VARIAC can control the voltage from 0 to 130% of the nominal voltage (l15 V) and
has a maximum capacity of 7.5 A of current. Additional power input is available by using
the electric heater. A second VARIAC controls the power input produced by the electric
heater.
3.2.3ı Instrumentation
The concrete slab and the hydronic heater are fully instrumented. The temperature
measurement is done with 20 type T thermocouples. It is possible to observe the position
ofsome thermocouple during the construction process in Figure 3.2-3. A 2-D schematic
75ˇ
drawing shows the thermocouple position in Figure 3.2-7. The coordinates where the
thermocouples are placed in reference to the central pipe are shown on Appendix C.
H4 H3 H2
• 9
.1
-11
·12 .6
0 18 19 20 0 0
013 -8
• 14
SECTION A - A'
Figure 3.2-7: Position ofthermocouple at the concrete slab
Only thermocouples 4 and 3 (inlet and outlet ofthe embedded piping) are
calibrated.. The thermocouples were calibrated to ± O.OsoC. The calibration process can
be found in appendix A. Unfortunately, there was no documentation about the calibration
for the rest ofthe thermocouples used in the concrete slab. Therefore, the expected
accuracy for them is ± O.S°c.
76
The power used. by the circulation pump and the heater is also measured with watt
transducers installed in line with each component. Originally, the flow rate measurement
was intended only as a control parameter only and was not important for the power
calculation. BU4 due to uncontrolled heat losses from the piping outside the slab, it is not
possible to use the measured value directly. Still, the power measured by the watt
transducers is used as a check. The heat transfer to the slab is estimated using the flow
rate reading available from the in-line rotameter and the inlet and outlet thennocouples.
The in-flow rotameter is a VTP-4020. This model provides a measurement with
accuracy of ±8% full scale (±O.8 Lpm or ±0.20 GPM). This rotameter is viscous
independent for fluids in the range of 1 to lOOcSt. The variation for the viscosity for the
Ethylene Glycol is between 1.85-4.55 cSt for these application. The rotameter cannot be
accessed during the experimentation process, so the flow rate is measured at the
beginning and end of the experiment.
One FLUKE Hydra II data logger is used to store the temperatures and reference
power. The time interval used was 5-minutes. The infonnation is downloaded and
processed after the test.
3.3 Snow making I Environmental Chamber
The environmental chamber was designed to study the snow melting phenomena
in interaction with the concrete slab. The system facilitates the study in the laboratory of
similar process occurring on a bridge deck. Since it is not possible to get fresh snow all
year round (it is even difficult to obtain during the winter season in Oklahoma), it is
necessary to have equipment to produce artificial snow when required. Once the snow is
77
obtained, it is also necessary for the study purposes to maintain and control the
environmental condition. This equipment can accomplish both tasks. Additionally, the
environmental chamber has a specific compartment in which the concrete slab is placed.
The main criteria of design for the environmental chamber were the required
dimensions. The first requirement was the ability to hold the concrete slab inside the
chamber. For this reason the working area ofthe snow chamber top and central section is
O.91m x O.91m (3ft x 3ft). The second requirement was that the size ofthe chamber
should fit in the available room. This limits the hei~t to 2.6Om (8.53ft).
The environmental chamber was constructed in three separate sections in order to
facilitate the movement to a different location ifnecessary. Each section has a specific
objective. The bottom section holds the concrete slab. The central section has the
equipment necessary to generate snow. The top section provides the necessary space for
the snow to be formed. The top and central sections have the same working area in order
to allow the formed snow to fall directly over the concrete slab. The bottom section has a
working area of 1.52m (5ft) by 1.52m (5ft) in order to allow the whole concrete slab and
the circulation equipment to fit into the chamber.
The front side ofthe environmental chamber has two access doors. The first is
located in the bottom section and is used to allow the concrete slab to be rolled into the
chamber. A door located in the central zone allows the operator to obtain digital pictures
from the inside ofthe chamber during the melting process. Further discussion will he
made in section 3.3.2. Both doors are insulated and closed during operation. The three
sections are shown in Figure 3.3-8.
78
Figure 3.3-8: Environmental chamber sections
3.3.1ı Top section
The top section is a l.oOm (3.25ft) extension of the central section. The
dimensions of the inside box are the same as in the middle section 0.9lm x 0.9lm (3ft x
3ft). In the roof ofthe top section is removable and contains a fluorescent light. The walls
were built with 3.8cm x 14.0cm (1.5in x 5.5in) pine wood construction on OAm (l6in)
centers and filled with l3ern (5in) thick extruded polystyrene insulation. The inside walls
were covered with 1.6 mm (0.0625 in) ofhard plastic sheets that were sealed with
silicone caulking to the surface ofthe insulation. This was done in order to provide an
79ˇ
jacket that prevent damage to the insulation.
The roof is made of extruded polystyrene since it has no associated load. The
inside of the roof is covered with a bard plastic sheet to prevent water infiltration into the
fluorescent light. The light was added to aid in photographic documentation of the
experiments. The light is a standard 1.2m (4ft) fluorescent shop light. A sheet of
Plexiglas was placed over the rectangular slot from the inside to reduce heat gain to the
chamber and let the light shine inside. The light is controlled to automatically tum on
during the picture acquisition only. This is discussed in section 3.6.
A O.127m (5in) diameter hole was made on the backside ofthe top section. It
holds a O.61m (2ft) long piece of polyvinyl chloride (PVC) pipe used as the air return
duct for the mechanical refrigeration system. Ifno mechanical refrigeration is needed this
hole can be blocked with PVC end caps. The mechanical refrigeration system will be
discussed later on in this section.
3.3.2ı Central section
The central section has the same cross-sectional dimensions as the top section.
The central section height is 1.6Om (63in). Together the central and top section provides
the room necessary for the snow to form. Longwill, et a1. (1999) found that the snow
making process require a height of2..lm (7ft) for the snow to be produced. Further
development and experimentation ofHockersmith (2001) shows that a better quality of
snow can be achieved with longer residence time for the water droplets. Consequently,
the total dimension ofthe tower is 2.6Om (8.53ft).
80
The central section contains the snow making equipment. Two holes were made
in the bottom of the sidewalls. On each side, a water nozzle uses the first hole and a
liquid nitrogen nozzle uses the second. The holes were filled with hard plastic pipes and
sealed with silicone. This prevents water from entering into the insulation or wood
damaging it. The hard plastic pipes also allow an easy introduction of the nozzles into
position. The snow making equipment will be discussed later on in this section.
The central section also contains the O.61m by 1.22m (2 ft by 4 ft) access window,.
which has two purposes. The first one is to allow the Qperator access to the chamber. The
window is hinged to allow access. The second purpose is to provide a clear view ofthe
inside of the chamber while the system is running. A smalll opening of 12.7em x 12.7em
(5in x 5in) was made in the clear cover. Thru this opening the digital camera's lens is
placed. This device is used to acquire digital pictures ofthe surface ofthe slab during the
snow-melting process.
During snow production the window is not in use, so a 12.7cm (Sin) extruded
polystyrene plug is used to cover it. On the other hand, during the melting process digital
pictures are required. So, an alternate plug that holds the space for a Kodak digital
camera replaces the standard plug. In order to prevent heat losses from the chamber the
lens opening is sealed with Polyurethane foam. The lens opening is shown in Figure
3.3-9.
81ı
Figure 3.3-9: aJ Lens opening. h) Digitalcamera space in alternateplug
3.3.3 Bottom section
The dimensions ofthe bottom section were established based on the requirements
for the concrete slab and the diffuser for the mechanical cooling system., 'fhis s~~on
includes two doors ofO.9lm (3ft) each in the front side to provide enougb room for
handling the concrete slab while loading. The material used for the doors areex:truded
polystyrene with a hard plastic cover for the inside, and sheet metal on the outside. The
rest of the body construction is similar to the other sections using 3.8cm. x 14.0cm (1.5in
x 5.5in) pine wood construction on O.3Om (1ft) centers and 13cm (5in) of extruded
polystyrene insulation. On th.e inside, the walls are covered with hard plastic and sealed
with silicone caulking. The floor ofthe section is reinforced with 3.8cm x 14.0cm (1.5in
x 5.5in) pine wood construction placed on 20.3cm (Sin) centers in order to support the
weight ofthe concrete slab.
A water collector and a drain hole were placed in the backside ofthe floor. The
collector is used to drain the melt water. The drain channel, which runs from the side to
side at the back ofthe chamber, is 2.54cm (lin) depth, 12.7cm (5in) wide. To remove the
82ı
The mechanical cooling diffuser is connected inside the bottom section by a
12.7cm (Sin) hole. A O.6lcm (2ft) long piece ofPVC pipe is inserted into the hole and
sealed with silicon caulking. The diffuser has a U shape and is connected in the central
section to the inlet pipe. The design of the diffuser is such that the velocity of the flow is
kept to a minimum in order to prevent forced convection effects. The diameter of the
diffuser is 12.7cm (5in) with 2..54cm. (lin) holes. If the mechanical cooling is not used or
the conditions don't allow the use ofthe diffuser, it can be removed and a PVC cap seals
the inlet.
Figure 3.3-10: Mechanical cooling diffuser
3.4 M.echanical cooling
The environmental chamber is provided with a mechanical cooling system. This
system is used to initially cool the environmental chamber to temperatures near freezing
conditions. It also recreates the environmental conditions ofthe slab after the snow has
been produced. The mechanical cooling apparatus is connected to the environmental
83ı
located in the top section. The elements used to produce the mechanical cooling are
discussed next.
The system works with two recirculating chillers. The units consist ofa Non-CFC
air-cooled refrigeration system, circulation pump, work area bath (7L), and a
microprocessor temperature controller. ThennoNeslab is the manufacturer of each unit
(www.thennoneslab.com). The first unit has a cooling capacity of 750W (@ O°C) (RTE740)
and the second unit has a cooling capacity of500W (RTE-140) (@ O°C). The
pumping capacity ofboth chillers is 151 Iprn @ 0 m head and 0 lpm @ 4.9 rn head. The
working fluid used is Ethylene Glycol with a volumetric concentration of 50%. Each of
the chillers is connected in parallel arrangement to two fin-tube heat exchangers. The
heat exchangers are located in an insulated box. It is made by 12.7cm (5 in) of extruded
polystyrene insulation. The outside is covered with 6mm (0.25 in) plywood. Inside the
box is covered with 1.6mm (0.0625 in) ofhard plastic.
The heat exchanger box is connected to the inlet ofthe environmental chamber
and to a fan box. The fan box is constructed in the same way as the heat exchanger box. It
is covered with 0.8mm (0.03 I 3in) sheet metal jacket. This fan box is used to supply air to
the environmental chamber. Consequently, the fan inside the box can produce 0.08 rnJ/s.
The shaft ofthe fan was mounted thru the wall. So, the motor is installed outside the fan
box in order to reduce as possible the heat gain by the system. The AC motor is a
120V/60Hz with power of 125 W.
84ı
connected to the top section ofthe environmental chamber. This duct was constructed
with 10.16cm (4 in) PVC pipe and insulated in the outside with R-25 fiberglass bat
insulation. The overal1length ofthe duct is 4.6m (15 ft). A schematic drawing ofthe
system is shown in Figure 3.4-11.
I
Air return duct
Snow MakJng
Chamber
'- "'e ICl.l::! I Jltro~er. Tanl<
To Air Source .J1----=::::=r......r11Iı _ ~'- To Ir Source
To Water Sourceı
Concrete Slab Space
Diffuser
Heat
):::: ::::( Exchangers
I~==~"'~~~
L..--+
...I====----=-=~
Air supply duct Chillers
Figure 3.4-11: Snow making machine schemillics
3.5 Snow making equipment
The snow making equipment was designed in order to recreates the natural
formation of snow. According to Hockersmith (2002), the natural snow is formed by
85ı
At lab scale, it is difficult to facilitate the formation of snow in a comparable short
period of time. To maximize this effect the water droplets are jetted to the top of the
chamber. This duplicates the traveled space for the particle. The traveled distance in the
chamber is 5.2m (16ft). To form the jets several techniques have been suggested.
Longwill, et al. (1999) proposed the use of a bottom nozzle array that mix air and water
streams inside the piping system. In this way it boosts the water droplets using the 3
nozzles ofllie array. Hockersmith (2002) showed that the operation of the system is
difficult and instead used two modified engine degreaser nozzles. In this approach water
is pressurized to the nozzle and controlled with a trim valve on the degreaser. This system
was very effective, but the trim valve was easily damaged. Longwill, et a1. (1999) also
showed that the temperature ofwater used for snow production has huge effect in the
snow quality. For this reason, water near freezing is used in the system.
The system originally used was modified to avoid the use ofthe degreaser trim
valve and to allow the use of pre-cooled water. The modification consists in connecting
directly the air and water inlets to the nozzle. In this way the use of the trim valve was
eliminated. Instead a globe valve was installed outside the nozzle body and directly on
the water inlet pipe of each degreaser. The assembly is then connected directly to a precooled
water tank. The water inlet is controlled with the ball valves before each nozzle
and the air is controlled using the regulating valve located in the control panel ofthe
snow machine. See Figure 3.5-12a.
86ı
nitrogen into the chamber. The tubing is introduced by two 1.27em (0.5 in) holes 12.7cm
(5 in) above the water injection nozzles. The holes are fitted with 1.27cm (O.Sin) PVC
pipe and sealed with silicone caulking. The tubin.g is covered with 2.54 em. of foam wrap
insulation. On the inside ofthe chamber the tubing is bend to direct the liquid nitrogen
stream to water particles stream as shown in Figure 3.5-12. This rapid cooling helps on
the snowmaking performance.
Figure 3.5-12: a) MOdified nozzle schematics, b) Nozzles outside chamber, c) Nozzles inside the chamber
3.6ı Data acquisition controller
In order to acquire the digital images ofthe concrete slab surface a digital camera
is used. A fluorescent light at the top ofthe chamber is switched on 2 seconds prior to the
camera shutter opening and 2 seconds after the shutter is closed. Moreover, the digital
picture has to be synchronized with the data logger data acquisition. All this effort is
necessary to reduce the uncertainty produced by the introduction ofunnecessary heat
from the light inside the chamber.
87ˇ
controller is timed to match the camera clock and reduce the total "on" time of the light
to 7 seconds in each 10 minutes interval. Further discussion on the controller system can
be found in Appendix B.
3.7 Concrete thermal conduction coefficient measurement apparatus
Two of the concrete samples made during the concrete slab construction were
used to calculate the conductivity of the concrete used. The samples were cut to produce
a 7.62 cm (3 in) diameter and 7.62 em (3 in) tall cylinders. After polishing the top and.
bottom sides each sample was independently tested using a Single-Sided Hot-Plate
apparatus similar to the one used in the ASTM CI 77-93 Standard.
This apparatus consists ofa heavy insulated box, an electric heater plate, and a
cooling plate. The insulation box is a 1.5 m x 1m x 1m box filled with extruded
polystyrene and works as primary and secondary guard for the hot plate. A perfectly
sized hole is made at the center ofthe insulation package were the sample, the heater
plate, and the cooler plate are introduced as a single set. The testing pack can be seen in
Figure 3.7-13.
88ı
Figure 3.7-13: CONcrete testiNg pack
The heater plate is located in the left side of the pack. The heater uses an electrical
resistance to produce the heat input. The sample and the plate are joined together by a
thin film of a highly conductive paste (with known conduC4mce) and held in place using
duct tape. The heater plate is also provided with two temperature sensors (thennistors) to
measure temperature symmetrically at the sample surface.
The cooling plate is located in the right side of the pack and it is used to keep the
temperature of the sample from overheating. The circulation fluid is filtered water that
runs in an internal channel of the plate. It is joint to the sample in the same way that the
heat plate. The cooling plate has also two temperature sensors (thermistor) to measure the
temperature symmetrically over th.e surface of the sample.
A 16-bit DAQ board connected to a PC computer does the data acquisition. The
data samples are taken every minute. Consequently, it is possible to obtain a detailed
temperature profile from each side ofthe sample and use it to calculate the conductivity
89ı
chapter.
3.8ı Specific heat calculation apparatus
One concrete specimen was prepared for experimentation. The sample had the
same characteristics as the one used for the conduction experiments. The only difference
is that this sample has been supplied with a thermocouple placed in the center ofthe
sample. The same thennocouple cable is used to place the sample in the correct position
while inside the calorimeter. This is done to avoid the introduction ofan alien object to
support the sample into the calorimeter (i.e. wire basket to hold the sample). See Figure
3.8-14.
Figure 3.8-14: Concrete sample/or Cp experimentation
The calorimeter is a metallic Dewar flask with capacity of2 liter . It is covered
with a 101.6mm (4 in) extruded polystyrene insulation. This is done in order to prevent
heat transfer from the environment. A blue foam cover is used at the top of the
90ˇ
sealed with caulking. See Figure 3.8-15.
Figure 3.8-1S: Calorimeter insulation
An important part of the experiment is the drying of. the sample. An electric oven
was used to dry the sample. The sample needs to dry for 4 hours at a rated tem.perature of
150°C (302 OF). After the oven-dried period, the sample must rest in a water bath for a
period of 48 hours. Finally, the sample is weighted and introduced into a boiling water
bath. The weighting process is done in order to make possible the measurement of th.e
amount of water carried by the sample during the heating in the water bath. After the
oven dry period, the wetting period, and the water bath period, the sample wa weighted
using a digital balance readable within ± O.lg. The digital balance is shown in Figure
3.8-16.
91ı
Figure 3.8-16: Digital balance
The temperature during the test period was measured at three different points. The
first measurement point was the water bath used to heat up the sample. The econd
measurement point was the concrete sample. The third measurement point was placed
inside the calorimeter to measure the temperature of the contained water. All the
temperatures were acquired with cooper-constantan (Type n calibrated thermocouples to
± o.osoe (± 0.09 OF). The temperatures were recorded using a HP Hydra II data
acquisition device.
92ı
EXPERIMENTATION AND RESULTS
4.1 Introduction
This experimental investigation is oriented to validate the two snow melting
models used to simulate and optimize snow-melting systems. These models are the Finite
Difference Rectangular Grid (FD-RG) and the Finite Volume Boundary Fitted Grid (FVBFG)
models.
Several tasks were required in order to use these models. First, a model of the
entire experimental apparatus (concrete slab and supporting frame) was developed. Since
the specific thermal properties of the concrete used in the slab were not known, additional
experimentation was necessary.. The parameters investigated were conductivity and
specific heat of the concrete. Third, modifications of the codes to accommodate the
environmental conditions inside the snow chamber and the snow characteristics within
the models were made. The final task involves the description of the initial conditions in
each model. Therefore, the way in which each model is initialized is discussed in this
chapter.
Two different scenarios were analyzed in order to validate the models: dry surface
and snow covered conditions. Since a dry surface condition is the simplest case, it was
93ı
model of the concrete slab. The second condition analyzed is the more complex scenario
and used to validate the mathematical model of the snow-melting process. The
experimental procedure used and the results obtained in each analyzed case are described
in this chapter.
4.2ı Thermal property measurements
Even with a correct implementation ofthe weather infonnation the model will
return wrong results without proper physical data for the materials used in the slab..
Transient analysis requires the identification of specific heat for the material used. In this
case the materials used are extruded polystyrene, plywood board, and concrete. The
specific heat and thermal conductivity of extruded polystyrene and plywood have much
lower volumetric specific heats than concrete and may be well characterized from the
literature. However, according to Whiting, et at (1978), concrete has a very wide range
of variation (specific heat: 600-1500 J/(kg K) and thermal conductivity: 0.8 to 1.8 W/(m
K» depending on the composition (density) and water content of the concrete used.
Therefore, it was necessary to measure the density, specific heat, and thermal
conductivity ofthe concrete used.
4.2.1ı Density measurement
The density of the concrete was established using the following correlation:
m
p=-ı (4.1)
V
Where:
94ˇ
p= Density, (kg/m3 or Ib/ft3
)
m = Mass~ (kg or lb)
A concrete specimen was prepared using one of the concrete samples described
shown in Figure 3.2-3. The dimensions ofth,e sample were measured dividing the sample
in quadrants. Therefore, two sets ofdiametric measurements on each side of the sample
(top and bottom) and four length measurements were acquired. The instrument used to
perform this task was a caliper readable within ±O.Olmm. The values obtained after the
measurements are shown in Table 4.2-1.
Quadrant Core Length (mm) Location Diameter (mm)
1 76.32 U(l-3) 76.04
2 75.84 U(2-4) 76.04
3 76.30 L(l-3) 76.02
4 76.83 L(2-4) 76.02
AVG 76.3 76.0
Table 4.2-1: Geometric measurement values/or concrete sample
The volume ofthe sample is therefore determined with the following equation
v= ~D2 L (4.2)
4
Where:
95ı
tr=: Pi constant, [dimension1ess]
D =: Diameter, [rom] or [in]
L =: Length, [rom] or [in]
The volume of the sample is 346572 ± I 17 rnm3
.
The concrete sample was oven dried during a period of4 hours. After this period,
the sample was weighed using a digital balance readable within ± O.lg. The oven dried
value obtained was: 725.4 ± O.lg and the saturated value:is 780.1 ± O.lg.
Using equation (4.1), the density value obtained is: 2093.3 ± 0.6 (kg/m3
). The
sample was retested under water saturation conditions and the obtained value is: 2250.8 ±
0.6 (kg/m3
).
4.2.2ı Specific heat of concrete
The test procedure used to calculate the specific heat value for the concrete used
in the slab followed the experimental procedure used by Whiting, et al. (I 978). This is the
same process used by the U.S. Anny Corps of Engineers Specification CRD-C-124-73
(Method ofTest for Specific Heat of Aggregates, Concrete, and Other Materials). The
procedure consists ofoven-drying the sample and then immersing it in a water bath at
ambient temperature for at least 48 hours. After this period the sample has reached
saturation and is weighed after drying the surface from remaiIring water. The sample is
then placed in a boiling water temperature bath (99.1 ± 0.5 °C) at the local barometric
pressure. The sample remains in the bath for at least 20 minutes. Then the sample and the
96ˇ
o - Carry over water.ı
S - Sample.ı
e - Water equivalent of calorimeter.ı
Whiting, et al. (1978) used a standard material (sapphire) as calibration material.
The next equation was used for this purpose:
(4.4)
Where:
M = Mass, (kg or lb).
Cp = Specific heat, [J/(kg K) or BTU/(lb OR)].
T = Temperature raise or fall, (K or OR).
Sub indexes:
e - Equivalent to water.ı
s- Sapphireı
b-Basket.ı
0- Carry over water.ı
1 - In calorimeter.ı
98
In the same way, the water equivalent ofthe calorimeter can be calculated by
using pure water as calibration material instead ofthe recommended Sapphire. Thus, a
simplified equation is required. This equation does not take into consideration the basket
required for the sapphire or the water carried over. Consequently, the water equivalent of
the calorimeter can be found using the next equation:
(4.5)
Where:
Me = Mass equivalent of calorimeter, (kg or Ib)
Me = Mass of calibration water, (kg or Ib)
Ml = Mass ofwater in calorimeter, (kg or Ib)
Tc = Temperature rise of fall of calibration water, (K. or OR)
T, = Temperature rise or fall ofwater in calorimeter, (K or OR)
For this experiment, the mass ofcalibration water used was 0.25 ±0.005 kg and
the mass ofwater in the calorimeter was 1.75 kg (3.861b). The temperature fall of the
calibration water was 38.32 ± 0.05°C and the temperature rise ofthe water in the
calorimeter were 5.17 ±0.05 °e. Therefore, using equation (4.5), the mass equivalent of
the calorimeter was found to be: 0.103 ± 0.060 kg (0.227. ± 0.13231b).
99ı
The mass ofthe carry-over water is calculated by measuring the sample before
and after the boiling period. The carry-over water was measured to be: 0.04 ± 2.'oE-4 kg.
Using equation 4.3, it is now possible to calculate the value ofthe specific heat
for the' concrete under SSD conditions. The obtained value is: 1036 ± 42 J/(kg K) [0.21 ±
0.01 BTU/(lb OF)].
According to Whiting. et aI. (1978) it is possible to establish a relationship for
specific heat ofthe concrete as function of specific heat under SSD conditions and the
actual moisture content. The relationship takes the form:
CPSSD +Cwr(Y-l) C = ---'===----:.;...~-~ (4.6)
p l+r(y-l)
W"sSD-WOD r= (4.7)
~SD
Where:
Cp = Specific heat ofconcrete at any moisture content, [J/{kg K) or BTU/(lb OR)]
CpSSD = Specific heat of saturated-surface dry (SSD) concrete, [J/(kg K) or
BTU/(lb OR)]
Cw = Specific heat ofwater, [J/(kg K) or BTU/(lb OR)]
y= SSD moisture content, (dimensionless)
100ı
(dimensionless)
WSSD = SSD weight ofconcrete specimen, (kg or lb)
WOD =Oven-dry weight of concrete specimen, (kg or lb)
Equation 4.6 can be used to estimate the new specific heat for the wet concrete
during a snow test. The value of gamma ty) corresponding to SSD moisture content was
calculated as: 0.0706 (7.06%). Then the specific heat will range between: 799 J/(kg K)
[0.191 BTU/(lb OR)] for oven-dry conditions to 1036 J/(kg K) [0.247 BTU/(lb OR)] for
saturated (SSD) conditions.
4.2.3 Thermal conductivity
This parameter is measured using Hot-Plate Apparatus similar to that described in
the ASTM C177-97 standard (ASTM, 2002). A detailed description ofthe apparatus used
can be found in Section 3.7.
Sample Thickness
N TION
llo llo " '" 6 . ,.. 6 II <l
a l>
·llo
a a II
GI I> ~ :i! .'" ll. ' f i=
VI <l <l VI a. ~ ~ A A
"
:::E: " A i I C:l l> c..>
9
" " a.
A
II
+
6
N +
Paste thickness 12 P.lste thickness. 2
Figure 4.2-1: Schematic ofheat conduction on concrete sample
101
calculated from the steady state model as follows:
(4.8)
Q" kp
Where:
ks = Concrete sample conductivity, [W/(m K) or BTU/(hr ft OF)].
Kp = Paste conductivity, [W/(m K) or BTU/(hr ft OF) ].
t1T =Differential temperature between cold and hot side, eC or OF).
Q" = Heat flux on hot side, [W/m2 or BTU/(hr :tt2)].
Ls = Thickness of concrete sample, (m or ft).
Lp = Thickness of paste sample, (m or ft).
The thennal conductivity for the paste is given by the manufacturer (0.416 ±
0.0042 W/m-OC). By using this value and measuring the dimensions, all the parameters
ofthe equation 4.8 are set. Two test samples were analyzed. The conductivity value
found was: 1.15 ± 0.0063 W/(m °C) [0.664 ±0.0036 BTU/(hr ft OF)] at oven-dry
conditions and 1.61 ± 0.0063 W/(m °C) [0.664 ± 0.0036 BTU/(hr ft OF)] at saturation
conditions following Rhodes (1978) criterion.
102
found in appendix D. A calibration done with a material (Stainless Steel 308) with known
thennal conductivity gave errors ofless than ±2% (Smith, 2003).
4.2.4 Water diffusion in concrete experiment
An experiment was perfonned in order to investigate the rate ofwater diffusion
into the concrete slab. One ofthe cylindrical concrete slab samples (Figure 3.2-3) was
used. The sample was inserted into a close-fitting Mylar tube. Then the top surface ofthe
concrete was sealed around the outer edge. so that water could be ponded on the upper
surface ofthe concrete without leaking around the edges. Then 160 mt ofwater was
placed on the upper surface ofthe concrete shib sample. A picture is taken every 5
minutes to observe the water diffusion process into the concrete sample.
After the water was placed on the sample, it was noticed that the porosity was
higher near the top ofthe sample than at the bottom. Some water quickly passed throu.gh
the upper surface to the space between the Mylar tube and the concrete sample as shown
in Figure 4.2-2. The percentage of the outer surface of the sample that was saturated was
estimated visually at 5-minute intervals. The estimates and an equation are plotted in
Figure 4.2-2.
103
90
80
70
60
I.,. 50
40
30
20
10 I--------jy• O.ooo2x" . 0.0395x' + 3.3087. ·0.35111--__--1
~·0.9996
0
0 to 20 30 40 50 60 70 80 90 100
tlme(mln)
I • 'Yo waler Poly. (% water) I
Figure 4.2-2: Water dijfllsion in concrete slab sa/llple
The fonowing equation is used to correlate the water content with time:
(4-9)
Where:
we = water content (%)
t = time, [min]
CII = 2xlO-4, [lImin3
]
C/2 = -3.96xlO-2, [lImin2
]
C13 =3.3087, [l/min]
104
Since the concrete sample and the concrete slab have the same thickness (l5Omm
or 6 in) the penetration depth at any time may be approximated as:
PD= we *Th (4-10)
100 s
Where:
PD = Penetration depth
Ths = Concrete slab thickness (0.152 m or 0.5 ft)
we = water content (non-dimensional)
The penetration depth is used to change the specific heat and thermal conductivity
values from dry to saturated conditions. This procedure is done on a cell-by-cell basis to
approximate the water absorption ofthec..oncrete during the snow-melting phase ofthe
experimentation.
Once the snow melting is finished the specific heat and thermal conduction of
saturated cells are returned to dry conditions. The process is assumed to occur on a cellby-
cell basis starting from the deepest row to the surface. The application ofthis
experiment is shown in section 4.5.3.
4.3 Mathematical model of the concrete slab
The mathematical model of the concrete slab is based on a center section as
shown in Figure 4.3-3. Two planes of symmetry (through the tube and between the tubes)
105
- - -------------------------------------------- ---------------------------------
------ ----
-------------------------------------------- ---------------------- -----------
-----------
that the fluid temperature remains approximately constant along the tube, and therefore a
2-D model is appropriat.e.
4.3.1 Heat conduction
During construction, the central section of the testing slab was reinforced with
wood, and the void space filled with extruded polystyrene. Therefore, the central section
became a very complex structure as shown in Figure 4.3-3a.
.. .. ., '. . I -. .
~ ~
'.. '. . . .. ;; ' · . . .' • •• ". 1" ".11 •• ". ". 4. 1· .... • .. .....
CONCRETE' .
.. '\ -._~~~:.R~E -',-1
• ... ... ..« .. . .•. ..... . .•. ....
.... .. .~ .. .. ': ... ... ..
~ · .' .. . . . .. .. ... ' . ... ..
WOOD -----------
EXTRUDED ~. - =-
IDEALIZED =:: POLYSTYRENE 8
o :: :: MATERIAL :: =
WOOD ----------Figure
4.3-3: Central Section. a) Construction model, b) Mathematieal representation
106
the combined properties ofwood and extruded polystyrene was created. The combined
heat conduction was calculated assuming an ideal 1-D heat transfer case. This is a
reasonable assumption since the horizontal wood layers have relatively low conductivity
and do not serve as effective fins.
The model can be represented as a resistance network as shown below:
Rw2
Figure 4.3-4: Resistance representation ofslab's base
Equation (4.11) represents the diagram in Figure 4.3-4:
R~q =Rw. +(_I_+_1_J-1 +R
W2 (4.11)
RSF RW2
Where:
Req = Equivalent resistance, [K/W or (OR hr)/BTU]
107
the combined properties of wood and extruded polystyrene was created. The combined
heat conduction was calculated assuming an ideal 1-D heat transfer case. This is a
reasonable asswnption since the horizontal wood layers have relatively low conductivity
and do not serve as effective fins.
The model can be represented as a resistance network as shown below:
Rw2
Figure 4.3-4: Resistllnce representation ofslab's base
Equation (4.11) represents the diagram in Figure 4.3-4:
Req =Rwl +(_I_+_1_J-1+Rw2 (4.11)
RSF RW2
Where:
Req = Equivalent resistance, [KIW or eR hr)/BTU]
107
RWl = Wood resistance on vertical section, [K/W or eR hr)/BTU]
RSF = Extruded polystyrene resistance, [K/W or (OR hr)/BTU]
The resistance value can be calculated using the following equation:
R=~ (4.12)
kA.
Where:
R = resistance, [K/W or eR hr)/BTU] .
L = Length, (m or ft)
k = Heat conduction coefficient, [W/(m K) or BTU/(ft OF)]
A = Cross-sectional area, (m2 or ft2)
The data for the different resistances is tabulated in Table 4.3-2:
RWI RsF RW2
LenJnh (m) 0.0191 0.2477 0.2477
Area (m2
) 0.1277 0.1103 0.0174
Thennal conductivity, rW/(m K)l 0.17 0.03 0.17
Heat conduction resistance, (KIW) 0.8798 74.8565 83.7390
Table 4.3-2: Data/or resistll1lces ill slab
The equivalent value for the heat conduction coefficient is calculated by using
equation 4.12:
108
Where:
keq = Equivalent conduction coefficient, [W/(m K) or BTU/(hr ft OR)]
Req = Equivalent heat resistance coefficient [KIW or (OR hr)/BTU]
L = Total conduction length, (m or ft)
A complete geometrical representation of the central section of the concrete slab
can be found in Appendix C. The heat conduction coefficient for the ideal composite
material is: 0.08 W/(m K).
4.3.2 Specific heat and density
An equivalent specific heat for the ideal material ofthe slab base can be found
using the concept of volumetric specific heat. Thus, the sum ofthe volumetric specific
heats of the different materials has to be equivalent to a volumetric specific heat ofthe
idealized material. The equation used is shown below:
(4.14)
Where:
p = Density, (kg/m3 or Ib/ft3
)
Cp = Specific heat, [J/(kg K) or BTU/(lb OR)]
109
eq = Equivalent for section
Wi = Wooden section top and bottom
W2 = Wooden vertical section
SF = Extruded polystyrene section
T = Total volume.
The values used for this calculation are resumed in Table 4.3-3:
Rwl RSF
,I
, Rw2
Density, (kg/m3
) 415.00 43.20 415.00
Specific heat, [J/(kg K)]
Volume, (m3
)
2385.0
0.0049
1220.0
0.0273
2385.0
0.0043
Total value, (JIK) 4849.89 1438.82 4256.03
Table 4.3-3: Data uSlldfor volumetric specific heat equivalent
The value found for the combined density and capacitance was: 288897 JIK. The
value ofdensity was found using the following relation:
(4.15)
'Where:
p = Density, (kg/m3 or Ib/ft3
)
110
eq = Equivalent for section
W = Wooden section top and bottom
SF = Extruded polystyrene