NUMERICAL AND EXPERIMENTAL 

STUDY OF IN-DRIFT HEAT AND MASS 

TRANSFER PROCESSES

Prepared for 

U.S. Nuclear Regulatory Commission 

Contract NRC–02–07–006 

Prepared by 

Kaushik Das1

Chandrika Manepally1

Randall Fedors2 

Debashis Basu1

1

Center for Nuclear Waste Regulatory Analyses 

San Antonio, Texas 

2

U.S. Nuclear Regulatory Commission 

Washington DC 

September 2011
ii 

ABSTRACT 

Water vapor flow and thermal environment near the emplaced waste can significantly affect the 

long-term repository performance. An interplay of chemical, thermal, hydrological, and 

geomechanical processes creates a complex environment that affects the water chemistry and 

water distribution pattern. Water source and distribution can affect the integrity of the 

engineered barrier system. Computational Fluid Dynamics (CFD) is an established numerical 

analysis method that has been used successfully to simulate ground water and vapor flow and 

thermal fields for a variety of applications. It provides an effective tool for quantitative analysis 

of the multiphase heat and mass transfer that takes place inside the drift. This report describes 

CFD tool development and validation activities that have been carried out at the Center for 

Nuclear Waste Regulatory Analyses. The general purpose CFD code ANSYS-FLUENT®

Version 12.1 (ANSYS, Inc., 2009) was used to simulate the multimode heat and mass transfer 

within the drift. Modules were developed to simulate water evaporation and condensation, as 

the standard CFD solver does not have a suitable model to represent these processes. These 

modules were developed as customized functions and were linked with the main solver as 

User Defined Function. This report describes the development of the customized functions and 

in particular, assessment of these add-on functions to evaluate their suitability in simulating the 

in-drift heat and mass transfer processes. Two different experimental studies were conducted 

for benchmarking the simulated data. The first experiment was carried out in a closely 

controlled environment, where a single condensation-evaporation cycle was generated within a 

rectangular space. Temperatures at specific locations and condensation rates were measured. 

Subsequently, a model was developed using ANSYS-FLUENT to replicate the experimental 

setup and the customized functions were used to calculate condensation and evaporation rates. 

Computed results were compared with experimental data to gain confidence in the developed 

customized modules. The second experiment used for validation was a 20-percent drift-scale 

model of the repository using electrically heated cylinders to create the effect of heated waste 

packages and dripping water to replicate seepage. The experimental investigation studied both 

uniform and nonuniform heat load distribution and its effect on two-phase water movement 

within the drift. Numerical model results indicate the presence of strong vertical convection 

above the waste packages and within the gap between the cylindrical packages. Computed 

and experimental temperatures across the domain were compared and exhibited similar trends. 

Based on the benchmarking exercise, the developed modules were found to be suitable for their 

intended use for in-drift heat and mass transfer simulations. 

Reference: 

ANSYS, Inc. “ANSYS-FLUENT® Version 12.1 User’s Guide.” Canonsburg, Pennsylvania: 

ANSYS, Inc. 2009. 
iii 

CONTENTS 

Section Page 

ABSTRACT ................................................................................................................................... ii 

FIGURES ..................................................................................................................................... vi 

TABLES ...................................................................................................................................... viii 

ACKNOWLEDGMENTS .............................................................................................................. ix 

1 INTRODUCTION ............................................................................................................... 1-1 

 1.1 Background ........................................................................................................... 1-1

 1.1.1 Relevant Physical Processes ................................................................... 1-1 

 1.1.2 Discussion on Previous Studies ............................................................... 1-3 

 1.2 Organization of the Report .................................................................................... 1-5 

2 NUMERICAL METHOD .................................................................................................... 2-1 

 2.1 Description of the General Purpose Solver ........................................................... 2-1 

 2.2 Technical Requirements of Customized Functions ............................................... 2-2 

 2.2.1 Nomenclature ........................................................................................... 2-2 

 2.2.2 Wall Condensation ................................................................................... 2-3 

 2.2.2.1 Source Term Specifications ..................................................... 2-6 

 2.2.2.2 Boundary Condition Specifications .......................................... 2-6 

 2.2.3 Wall Evaporation ...................................................................................... 2-7 

 2.2.4 Volumetric Condensation ......................................................................... 2-9 

 2.2.4.1 Equilibrium Model ...................................................................... 2-9 

 2.2.4.2 Nonequilibrium Model .............................................................. 2-10 

 2.2.4.3 Source Term Calculation ......................................................... 2-11 

 2.3 General Numerics ............................................................................................... 2-11 

 2.4 Verification Calculation ........................................................................................ 2-11 

 2.4.1 Description of the Numerical Model ....................................................... 2-12 

 2.4.2 Test Results for Nonequilibrium Model .................................................. 2-15 

 2.4.3 Test Results for Equilibrium Model ........................................................ 2-16 

 2.5 Summary ............................................................................................................. 2-17

3 CONDENSATION CELL EXPERIMENT AND VALIDATION ............................................ 3-1 

 3.1 Description of Water Transport Experiments ........................................................ 3-1 

 3.1.1 Test Chamber .......................................................................................... 3-1 

 3.1.2 Observations and Test Measurements .................................................... 3-3 

 3.2 Numerical Modeling of the Moisture Transport Experiments ................................ 3-5 

 3.2.1 Numerical Method .................................................................................... 3-5 

 3.2.2 Modeling Results ...................................................................................... 3-6 

 3.2.3 Comparison Between and Experimental and Computed Data ................. 3-9 

 3.3 Summary ............................................................................................................. 3-11

4 20-PERCENT SCALE MODEL OF A DRIFT .................................................................... 4-1 

 4.1 The 20-Percent Drift-Scale Experiment ................................................................. 4-1 

 4.1.1 Design of the Experiment ......................................................................... 4-1 

 4.1.2 System Monitoring ................................................................................... 4-3 

 4.1.3 Experimental Conditions .......................................................................... 4-4 

 4.1.4 Experimental Results ............................................................................... 4-6 

 
iv 

CONTENTS (continued) 

Section Page 

4.2 The 20-Percent Drift-Scale Modeling ................................................................................ 4-6 

 4.2.1 Model Development ................................................................................. 4-6 

 4.2.2 Model Results ........................................................................................ 4-10 

 4.2.2.1 Simulation of Dry Cases .......................................................... 4-10 

 4.2.2.2 Simulation of Moist Cases ....................................................... 4-18 

 4.3 Summary ............................................................................................................. 4-24

5 SUMMARY ........................................................................................................................ 5-1 

6 REFERENCES .................................................................................................................. 6-1 
v 

FIGURES 

Figure Page 

2-1 Schematic of the Condensation Model .......................................................................... 2-4 

2-2 Schematic of the Evaporation Model ............................................................................. 2-8 

2-3 Schematic for Condensation of Humid Air Over Flat Plate .......................................... 2-12 

2-4 Air Mass Fraction Variation in the Flow Field ............................................................... 2-13 

2-5 Comparison of Analytical and Computed Results for Condensation Mass Flux for 

 Inlet Velocity = 0.1 m/s [0.33 ft/s] Using Nonequilibrium Model ................................... 2-13 

2-6 Comparison of Analytical and Computer Results for Condensation Mass Flux for 

 Inlet Velocity = 1 m/s [3.28 ft/s] Using Nonequilibrium Model ...................................... 2-13 

2-7 Relative Humidity Contours for Simulation Using Nonequilibrium Model and Inlet 

 Velocity = 1 m/s [3.28 ft/s] ............................................................................................ 2-14

2-8 Air Mass Fraction Variation in the Flow Field With Multiphase Flow and Inlet 

 Velocity = 1m/s [3.28 ft/s] ............................................................................................. 2-14

2-9 Comparison of Analytical and Computed Results for Condensation Mass Flux for 

 Inlet Velocity = 0.1 m/s [0.33 ft/s] Using Equilibrium Model ......................................... 2-14 

2-10 Comparison of Analytical and Computer Results for Condensation Mass Flux for 

 Inlet Velocity = 1 m/s [3.28 ft/s] Using Equilibrium Model ............................................ 2-15 

2-11 Relative Humidity Contours for Simulation With Multiphase Species Transport and 

 Inlet Velocity = 1 m/s [3.28 ft/s] .................................................................................... 2-15 

3-1 Side View Schematic of Laboratory Test Chamber ........................................................ 3-1 

3-2 Side View Photographaboratory Test Chamber ............................................................. 3-2 

3-3 Top and Side Views Schematic of Laboratory Test Chamber. All Dimensions Are In 

 Inches. T1 to T9 Are Thermocouples ............................................................................ 3-3 

3-4 Experimental Results Showing the Change in the Condensation Rate As a Function 

 of the Difference Between the Heated Water and Cold Plate Temperatures ................. 3-5 

3-5 Computational Grid and Boundary Condition ................................................................. 3-6 

3-6 Velocity (m/s) Contours Using Nonequilibrium Model (Without Volumetric 

 Condensation) [1 m/s = 3.28 ft/s] ................................................................................... 3-7 

3-7 Velocity (m/s) Contours Using Equilibrium Model (With Volumetric Condensation 

 [1 m/s = 3.28 ft/s] ............................................................................................................ 3-7 

3-8 Temperature (K) Contours Using Nonequilibrium Model (Without Volumetric 

 Condensation) [°F = 1.8 °K × −459.4] ............................................................................ 3-8 

3-9 Temperature (K) Contours Using Equilibrium Model (With Volumetric Condensation) 

 [°F = 1.8 °K × −459.4] .................................................................................................... 3-8 

3-10 Relative Humidity Contours Using Nonequilibrium Model (Without 

 Volumetric Condensation) .............................................................................................. 3-9 

3-11 Relative Humidity Contours Using Equilibrium Model (With 

 Volumetric Condensation) .............................................................................................. 3-9 

3-12 Interphase Mass Transfer (kg/m2

-s) [1 kg/m2-s = 0.205 lb/ft2-s] Due to Volumetric 

 Condensation Using Multiphase Flow .......................................................................... 3-10 

3-13 Comparison of Experimental and Values of Condensation Rate Using Nonequilibrium 

 Model (Without Volumetric Condensation) ................................................................... 3-10 

3-14 Comparison of Experimental and Computed Values of Condensation Rate Using 

 Equilibrium Model (With Volumetric Condensation) ..................................................... 3-11 

3-15 Comparison of Experimental and Computed Values of Temperature Using 

 Nonequilibrium Model (Without Volumetric Condensation) .......................................... 3-12 
vi 

FIGURES (continued) 

Figure Page 

3-16 Comparison of Experimental and Computed Values of Temperature Using 

 Equilibrium Model (With Volumetric Condensation) ..................................................... 3-12 

4-1 Photographs of the 20-Percent Scale Drift Natural Convection and Cold-Trap 

 Laboratory Model ........................................................................................................... 4-2 

4-2 Schematic Drawing of the 20-Percent Drift-Scale Laboratory Model. Longitudinal 

 Cross Section Shows the Length of the Tube and Waste Package and Their Location 

 Relative to the Hot End Wall .......................................................................................... 4-3 

4-3 Photograph of Moist Area on Top of Waste Package A During Experiments 

 With Moisture ................................................................................................................. 4-4 

4-4 (a) Schematic Showing the Location of Axial Planes Selected for Placing 

 Thermocouples and (b) Example Location of Thermocouples at a Cross Section ........ 4-5 

4-5 Observed Average Cross Sectional Air-Temperature Profile for All the Dry Cases 

 With Different Heat Load Distributions Between Individual Waste Packages ................ 4-7 

4-6 Schematic of the Three-Dimensional ANSYS-FLUENT® Model Showing 

 Relevant Dimensions ..................................................................................................... 4-8 

4-7 Schematic of the Three-Dimensional ANSYS-FLUENT® Model Showing Dimensions 

 and Boundary Conditions ............................................................................................... 4-8 

4-8 Volumetric Mesh for the Computing Domain Used in the ANSYS-FLUENT® Model ..... 4-9 

4-9 Highlight of Meshing Over Waste Packages and End Caps Used in the 

 ANSYS-FLUENT® Model ............................................................................................... 4-9 

4-10 Airflow Streamlines for a Uniform 50 W [170.6 BTU/hr] Heat Load ............................. 4-11 

4-11 Schematic of Cross Sectional Planes .......................................................................... 4-12 

4-12 Velocity Vectors and Axial Velocity Contours for the Dry Simulation With a 

 Thermal Load of 50 W .................................................................................................. 4-12 

4-13 Velocity Vectors and Axial Velocity Contours for the Dry Simulation With Thermal 

 Loads of 75-25-25-75 W .............................................................................................. 4-13 

4-14 Velocity Vectors and Axial Velocity Contours for the Dry Simulation With a Thermal 

 Load of 75 W ................................................................................................................ 4-13 

4-15 Velocity Vectors and Axial Velocity Contours for the Dry Simulation With Thermal 

 Loads of 80-60-40-20 W .............................................................................................. 4-14 

4-16 Temperature Contours Dry Simulations ....................................................................... 4-15 

4-17 Temperature Distribution on Waste Package Surfaces ............................................... 4-16 

4-18 Cross Section Averaged Temperature Distribution Along the Axial Distance .............. 4-17 

4-19 Temperature Distribution on Waste Package Surfaces ............................................... 4-18 

4-20 Comparison of Experimental and Computed Temperature Data for all Dry Cases 

 With Uniform and Nonuniform Heat Load Distributions................................................ 4-19 

4-21 Axial Location of the Relative Humidity Probes for the Moist Case Experiments ........ 4-22 

4-22 Relative Position of 12 Thermocouples and Relative Humidity Sensors in the 

 Middle Vertical Plane ................................................................................................... 4-20

4-23 Cross Section Averaged Temperature Distribution Along the Axial Distance 

 for Moist Cases ............................................................................................................ 4-21 

4-24 Comparison of Experimental and Computed Temperature Data for All Moist Cases 

 With Uniform and Nonuniform Heat Load Distributions................................................ 4-21 

4-25 Relative Humidity Contours Obtained From Moist Simulations ................................... 4-22 
vii 

FIGURES (continued) 

Figure Page 

4-26 Cross Section Averaged Relative Humidity Distribution Along the Axial Distance 

 for Moist Cases ............................................................................................................ 4-23 

4-27 Comparison of Experimental and Computed Relative Humidity for the Moist Cases 

 With Uniform 50 W [170.6 BTU/hr] and Nonuniform 75-25-25-75 W 

 [266-85-85-266 BTU/hr] Thermal Load Distributions ................................................... 4-24 
viii 

TABLES 

Table Page 

3-1 Test Results and Test Conditions. Temperature (Temp.) Sensor Locations 

 (T3–T9) Are Shown in Figure 3-3 ................................................................................... 3-4 

4-1 Thermal Loads for the 20-Percent Scale Experiments................................................... 4-6 

 
ix 

ACKNOWLEDGMENTS

This report describes work performed by the Center for Nuclear Waste Regulatory 

Analyses (CNWRA®) and its contractors for the U.S. Nuclear Regulatory Commission 

(USNRC) under Contract No. NRC–02–07–006. The activities reported here were performed 

on behalf of the USNRC, Office of Nuclear Material Safety and Safeguards, Division of 

High-Level Waste Repository Safety. This report is an independent product of CNWRA 

and does not necessarily reflect the view or regulatory position of USNRC. The USNRC staff 

views expressed herein are preliminary and do not constitute a final judgment or 

determination of the matters addressed or of the acceptability of any licensing action that may 

be under consideration at USNRC. 

The authors gratefully acknowledge T. Wilt for his technical review, E. Pearcy for his 

programmatic and editorial reviews, and B. Street for her administrative support. The authors 

would also like to acknowledge S. Svedeman, S. Green, and D. Bannon for their help in the 

design, fabrication, and performance of the test runs of the condensation cell. B. Bell and 

A. Walavalkar from ANSYS, Inc. are acknowledged for their advice and assistance in setting 

up the simulation cases and development of the User Defined Functions used in this study. 

M. Muller assisted with postprocessing results and developing graphics. 

QUALITY OF DATA, ANALYSES, AND CODE DEVELOPMENT

DATA: All CNWRA-generated original data contained in this report meet quality 

assurance requirements described in the Geosciences and Engineering Division Quality 

Assurance Manual. Sources of other data should be consulted for determining the level of 

quality of those data. The work presented in this report is documented in Scientific 

Notebook No. 841E (Das, 2011). 

ANALYSES AND CODES: The general purpose computational fluid dynamics simulation code 

ANSYS-FLUENT® Version 12.1 (ANSYS, Inc., 2009) was used to generate results for this 

report and is controlled in accordance with the CNWRA software quality assurance procedure 

Technical Operating Procedure–018, Development and Control of Scientific and Engineering 

Software. Postprocessing was done using TECPLOT 360 Version 2009 (TECPLOT, 2009). 

References:

ANSYS, Inc. “ANSYS-FLUENT® Version 12.1 User’s Guide.” Canonsburg, Pennsylvania: 

ANSYS, Inc. 2009. 

Das, K. “Unsaturated Zone Flow and Transport.” Scientific Notebook No. 841E. San Antonio, 

Texas: CNWRA. 2011. 

TECPLOT, Inc. “TECPLOT 360 Version 2009.” Bellevue, Washington: TECPLOT, Inc. 2009. 
1-1 

1 INTRODUCTION 

The in-drift repository environment is a result of a complex interaction between a number of 

thermal, hydrological and chemical processes. Numerical modeling techniques such as 

Computational Fluid Dynamics (CFD) that take into account the mutual interactions among 

these processes, can be useful to study such a multifaceted physical problem. The commercial 

off-the-shelf CFD package ANSYS-FLUENT® version 12.1 has been used for developing a 

baseline model for simulating the in-drift air flow and multimode heat transfer. In addition, 

customized User Defined Functions (UDFs) have been developed to include hydrological 

processes such as evaporation and condensation. This report summarizes the CFD model 

development, verification and benchmarking activities that have been carried out at the Center 

for Nuclear Waste Regulatory Analyses (CNWRA®). The verification and benchmarking 

calculation was supported by (i) a study of the condensation-evaporation cycle in a 

rectangular enclosure and (ii) a 20-percent drift-scale experiment of the proposed repository at 

Yucca Mountain. A brief description of these experiments is also provided in this report. 

1.1 Background 

This section presents a description of the in-drift physical processes that affect air flow and 

moisture movement. This section also provides a review of past experimental and modeling 

activities used in exploring the near-field in-drift thermohydrology. 

1.1.1 Relevant Physical Processes 

The water mass transfer and redistribution process within a drift is generally referred as 

cold-trap process. It is integrally linked with the thermally induced buoyancy driven natural 

convection of air within the drift as this is the primary mechanism of mass transfer. The 

cold-trap process is also fundamentally related to the interphase heat transfer processes at 

hot and cold surfaces, which act as the initial and terminal destination of water vapor within 

the system. 

The proposed emplacement of high-level waste in drifts will significantly elevate the 

temperatures of the drift environment. Convection cells caused by temperature gradients along 

drifts will lead to the movement of air and water vapor. In the cold-trap process, water 

evaporates at hotter locations, is carried in the vapor phase by convective air flow, and 

condenses at cooler locations. Elevated relative humidity combined with deliquescence may 

lead to liquid phase water contacting waste packages initially in peripheral zones of the 

repository where the drift-wall temperature may not exceed the boiling point. As time proceeds 

and the thermal perturbation decays, this zone with potential deliquescence migrates inward 

toward the center of the repository. The geometry of the components of the engineered barrier 

system (e.g., waste package and support, drip shield, and invert) and interaction with the 

wallrock complicates the problem of simulating air flow and condensation associated with the 

cold-trap process. 

Understanding convection in drifts is the first step for understanding the cold-trap process. Axial 

convection in the drift and latent-heat transfer will dampen axial temperature gradients. 

Offsetting this dampening is the heat flux from emplaced waste, which magnifies the 

temperature gradient between hot and cold locations. The intimate linkage of in-drift natural 

convection and condensation to heat and mass transfer in the host rock complicates modeling 
1-2 

efforts and requires the coupling of CFD and fluid flow in porous media models. The coupled 

simulations of porous media flow and in-drift flow using CFD codes will provide a detailed 

understanding of the involved physics. This coupling, however, is difficult because of the 

different time scales of the flow processes in the host rock and the drift. Porous media flow and 

heat transfer in the host rock has a much longer time scale, which ranges from months to years 

compared to the time scale of the convective flow in the drifts, which ranges from minutes to 

hours. Similarly, from Reynolds number (ratio of inertial to viscous forces in the flow space) 

considerations, independent CFD and porous media models are needed. Approximations of 

inertial flow effects at solid boundaries inherent in Darcy-based porous media models may not 

be valid for gas phase flow in large openings, such as the air spaces in drifts. These limitations 

have led researchers (Salari, 2005; Buscheck, et al., 2003) to develop independent models for 

the in-drift and surrounding porous host rock that are coupled through explicit exchange of 

boundary condition information between models at the interface. In general, separate CFD 

models are utilized for the large air spaces in drifts and porous media models for the fractured 

host rock. The question of explicit manual linkage or interactive dynamic linkage between these 

two approaches is a secondary issue at this point because of a larger question related to the 

sparseness of data to support any model results. 

Conductive, radiative, convective, and latent-heat transfer processes all influence temperature 

estimates in the emplacement drifts. For all the solid objects present, including those of the 

engineered barrier system, conduction plays an important role. Radiation heat transfer is 

significant in the air space between the drift wall and engineered barriers. Natural air 

convection above and below the drip shield will lead to cross sectional and axial air flow that 

enhances heat dissipation from the heated waste packages. Individually, conduction, natural 

convection, thermal radiation, and interphase latent heat-transfer processes are reasonably well 

understood. The combined effect of all these modes of heat transfer in a geometrically complex 

underground environment, however, is still considered an area of active research. A number of 

challenges exist in the numerical modeling of such a system, such as the proper boundary 

condition at the porous host rock surface, accommodating a complex geometry, proper 

modeling of all the heat transfer processes, etc. As mentioned previously, the most well 

accepted method of determining the boundary condition is to use estimated temperature 

gradients using analytical heat transfer or thermohydrological porous media models to provide 

the boundary conditions to the CFD model (Buscheck, et al., 2003; Danko and Bahrami, 2004). 

A number of simplifications in the geometrical configuration, based on their relative importance 

on the heat transfer process, are done to develop a CFD model that can provide useful results 

while using reasonable computing resources. With the intent to retain the most relevant 

physics, approximations and assumptions are also used to simplify the conjugate heat transfer 

model, which may lack generality in terms of the range of application but provides valuable 

insight of the in-drift conditions in the proposed repository. This report describes the 

development of a model that is intended for simulating the in-drift environment. This report also 

presents discussion of the assumptions and approximations made during the model 

development process regarding boundary condition, geometry and pertinent heat transfer. 

The general temperature and flow variation inside a drift are understood in terms of drift-scale 

{i.e., 1 km [0.6 mi]} and local-scale {i.e., <100 m [<328 ft]} variations. The local-scale variations 

are encountered in the immediate vicinity of a waste package and are generated mostly as a 

response to the magnitude of the thermal load of a single waste package or a series of 

neighboring waste packages. For example, local-scale variations are observed in the vertical 

natural convection pattern and the relative humidity around a waste package. Such local-scale 

convection cells can create a region of elevated or reduced temperature, reduced viscosity and 
1-3 

a preferential condensation zone within a drift. Drift-scale variations are broad patterns 

observed in the drift that are driven by the combined presence of a large number of waste 

packages within a drift. For example, the temperature difference between waste packages and 

drift wall will create a drift-scale convection that will tend to convect air in the axial direction. 

Similarly, heat load variation between two sets of waste packages that are placed far apart 

within a drift will start a convection current between the cooler and hotter region. It is clear that 

both local- and drift-scale variations contribute to the overall convection pattern and impacts 

moisture movement. A numerical model intended for resolving the in-drift flow and heat transfer 

processes should be capable of resolving both drift- and local-scale variations. 

The cold-trap process involves transport of water vapor from a warmer to a cooler location. 

Available liquid water on a warmer surface will evaporate, and the evaporated water vapor will 

convect with air. The water vapor can condense to form liquid water in two different ways. If 

the moist air comes in contact with a relatively cooler surface, water vapor will condense on the 

cold surface, forming liquid water and will release the latent heat of condensation to the surface. 

This process is called surface condensation. On the other hand, if the convecting moist air 

enters into a relatively cooler fluid zone, water vapor will condense to form mist. In this case, 

the latent heat will be released to the surrounding air water vapor mixture. This process is 

termed volumetric condensation. This report provides details of the technical activities related 

to formulating, implementing and validating the evaporation condensation model that was used 

to supplement the baseline ANSYS-FLUENT solver. 

1.1.2 Discussion on Previous Studies 

A number of studies to understand the in-drift heat and mass transfered have been done by the 

U.S. Department of Energy (DOE) and the U.S. Nuclear Regulatory Commission (NRC). This 

section provides a very brief overview of these activities. A detail review of the literature is 

provided by Manepally, et al. (2006). 

A comprehensive report by Bechtel SAIC Company, LLC (2004) described the DOE analysis 

pertaining to in-drift convection and condensation. DOE used two types of models to address 

in-drift convection and moisture redistribution: CFD and thermohydrological porous media 

models. CFD models focused on understanding in-drift processes and estimating the 

contribution of condensation to the flux of water entering the invert. Thermohydrological porous 

media models focused on incorporating the effect of in-drift processes on estimating 

temperature and relative humidity at the waste package and drip shield. To obtain the flow 

fields, DOE used the CFD code FLUENT®, which solves the time-averaged Navier-Stokes 

equations including turbulence and the energy equation including thermal radiation (Bechtel 

SAIC Company, LLC, 2004). Temperature variation along a drift in the surrounding host rock is 

included by incorporating 5 m [16 ft] of host rock (thermal model only) in the simulations. 

Two- and three-dimensional models are used to focus on issues at different scales, including 

representation of the drip shield, waste packages, and the invert in the grids. In the 

three-dimensional convection-only model, a 71-m [234-ft] portion of a drift (14 waste packages) 

is used as an analog for the much larger (~500-m [1640-ft]) emplacement half-drifts because of 

computational limitations. The two-dimensional representation uses the line-averaged waste 

package thermal loading and is used to gain confidence that radiation and turbulence in the 

gas-phase flow are accounted for in the three-dimensional model. A scaled down version of the 

repository illustrated the effect of including in-drift convection and moisture redistribution. A 

three-dimensional porous media model with three short drifts was simulated for three different 
1-4 

cases: (i) no axial transport in the drift or invert; (ii) a low value of a binary diffusion coefficient 

for axial transport, along with using pseudo-permeability and effective heat conductivity and 

adding a bulkhead; and (iii) a high value of a binary diffusion coefficient for axial transport, along 

with using pseudo-permeability and effective heat conductivity but no bulkhead (moisture and 

heat were allowed to escape to unheated portions of the drift). Bechtel SAIC Company, LLC 

(2004) concluded that the effect of convection and moisture redistribution on relative humidity, 

temperature, and rock saturation is encompassed in the broad range of uncertainty driven by 

percolation and host rock thermal conductivity. 

The effect of natural ventilation through access tunnels and shafts in removing moisture from 

the drift and drift wall is well recognized . The concept is that moisture is removed from the 

emplacement drifts by exchange of air with the ground surface. The air at the ground surface 

generally has a low relative humidity at ground surface temperatures, and consequently an even 

lower relative humidity at emplacement drift temperatures. Danko and Bahrami (2006, 2003), 

however, have simulated the effects of natural ventilation by including the air pressure 

differences between drifts and the ground surface to approximate the air exchange in their 

model. Additional modeling using integrated in-drift and host rock models discussed in 

Birkholzer and Danko (2005), which may more realistically track moisture redistribution in the 

reflux zone during the thermal peak and continue after the wallrock temperatures have dropped 

below boiling, may be used to further study the effect of natural ventilation on moisture removal. 

A number of independent studies performed at the CNWRA (Fedors, et al., 2004, 2003a,b) have 

evaluated the combined effect of in-drift natural convection with thermal radiation, conduction, 

and latent heat transfer using numerical (CFD) and laboratory models. These efforts have 

resulted in (i) development of increased understanding of the cold-trap process, (ii) identification 

of heat transfer processes to be explicitly represented in the CFD models, and (iii) experimental 

data not readily available in the literature that can be used to validate CFD models. 

A prototype benchtop laboratory experiment was developed, details of which are given in 

Fedors, et al. (2003a) to provide measurement data for supporting models of air flow and 

moisture redistribution in nonventilated heated tunnels. The initial modeling effort did not match 

the experimental results with a high level of accuracy (Fedors, et al., 2003a). Two important 

processes lacking in early modeling efforts were identified: (i) the CFD simulations of the 

benchtop experiment did not explicitly track moisture movement, including phase change and 

latent-heat transfer and (ii) the simulations did not adequately track heat transfer at solid 

boundaries. The first reason identified for the poor match between the model and the 

experiment (the lack of a phase change model in the CFD code) was fixed by developing a 

Moisture Transport Module for FLOW-3D (Flow Science, Inc., 2005) that accounts for the mass 

transfer and latent heat transfer associated with the evaporation and condensation of water in 

the drift (Green, et al., 2004). CFD results obtained using the Moisture Transport Module better 

represent the expected relative humidity near the heat source. The low measured values of 

relative humidity near the heat source invalidate the previous assumption of 100 percent relative 

humidity used prior to implementation of the Moisture Transport Module. A second reason 

suggested for difficulties in matching measured temperatures near the heat source in the 

benchtop cold-trap experiment was the inadequacy of standard heat transfer models at 

interfaces of solids and air (Fedors, et al., 2003a). Later testing of heat transfer models 

simulated using detailed grids within the boundary layers at the solid-gas boundaries helped 

reduce uncertainty caused by heat transfer models (Fedors, et al., 2004). 
1-5 

Additionally, two other experiments were developed to provide measured data to support the 

parameters used in the CFD models. The first involved a small condensation cell intended to 

provide data for moisture redistribution using a geometry that allows tight control of conditions. 

Tests were conducted to measure (i) water evaporation and water condensation rates in a 

natural convection flow and (ii) the amount of water transported from the water source to the 

condensation plate. The measured water transport rates were used to validate the 

condensation-evaporation model developed for simulating the cold-trap process within 

ANSYS-FLUENT. Details of the experiment and the comparison of experimental and numerical 

data are presented in Chapter 3. 

The second experiment used a 20-percent scale model of the emplacement drift with four 

analog waste packages providing a heat source in a long pipe. A detailed description of the 

20-percent drift-scale experiment and the related CFD model is provided in Chapter 4 of this 

report. The valuable insights learned from the prototype benchtop experiment guided many of 

the design features of the 20-percent scale experiment. Green and Manepally (2006) had 

developed customized routines for simulating the moisture transport and thermal radiation 

models that were incorporated into a commercial CFD package FLOW-3D (Version 9.0). 

Details of this modeling exercise are reported in Manepally, et al (2006). 

1.2 Organization of the Report 

This report describes activities related to development, testing and benchmarking of a CFD 

modeling approach used to simulate in-drift heat and mass transfer including the cold-trap 

process. The commercially available CFD package ANSYS-FLUENT Version 12.1 is used as 

the platform for model development. Customized UDFs are developed to incorporate the 

condensation-evaporation modeling capability within the baseline ANSYS-FLUENT solver. 

Chapter 2 describes the mathematical basis along with approximations and assumptions for the 

condensation-evaporation modeling equations. A two-dimensional validation calculation is also 

described in Chapter 2. Flow of moist air over a cold plate and resulting condensation is 

simulated using the baseline ANSYS-FLUENT solver and customized UDFs. Results presented 

include general flowfield information and comparison of computed results with a closed form 

analytical solution. Chapter 3 describes a combined experimental and computational study to 

investigate the evaporation condensation process in a rectangular enclosure. The experiment 

measured temperature at specific locations in the domain along with condensation rate. The 

numerical model simulated the experimental setup in a two-dimensional domain. Results show 

the velocity, temperature and relative humidity distribution. A comparison of experimental and 

simulated data for condensation rate and temperature is also presented. Chapter 4 describes 

the laboratory and numerical model results of the intermediate scale experiment (20-percent 

drift scale). Experiments and computations were carried out for uniform and nonuniform thermal 

loadings. Numerical results are presented for temperature, relative humidity and velocity 

distribution in the flow field. A comparison of measured and computed temperature and relative 

humidity data are also made. The summary of this report is provided in Chapter 5. 
2-1 

2 NUMERICAL METHOD 

Radioactive decay heat creates a natural circulation in the drift and heat dissipation takes place 

through conduction, convection and radiation. Additionally, the decay heat also initiates a 

condensation-evaporation cycle that involves latent heat exchange between phases. As 

mentioned in the previous chapter, simulation of in-drift thermal environment requires 

proper modeling of these coupled heat and mass transfer processes. Numerical models 

discussed in the report were developed using the general purpose commercial computational 

fluid dynamics (CFD) package ANSYS-FLUENT® Version 12.1 (ANSYS, Inc. 2009a,b). 

Though the baseline solver had adequate capabilities to model the heat transfer including 

radiation and mass transfer problem, development of customized functions was necessary to 

include the evaporation and condensation processes in the simulations. Capabilities of the 

baseline solver, technical basis for developing the customized functions, and their use are 

described in this chapter. 

2.1 Description of the General Purpose Solver 

ANSYS-FLUENT is a general purpose computational fluid dynamics (CFD) software package 

developed by ANSYS, Inc. It is widely used for design and analysis of applications ranging 

from aircraft components to sporting gear and garments. It is among the leading computer 

aided engineering tools that engineers use for flow, chemistry and heat transfer simulations 

in a continuum, though it also has some added modules for magnetohydrodynamics and 

discrete particle problems. To facilitate model development and pre and postprocessing 

activities, a number of additional software packages are also available. Three of these 

additional packages were used to build and analyze the models described in this report: 

(i) ANSYS-DesignModeler (ANSYS, Inc., 2009c), a computer aided design package for 

building geometries; (ii) ANSYS-Meshing (ANSYS, Inc., 2009d), used for mesh generation; 

and (iii) ANSYS-CFDPost (ANSYS, Inc., 2009e) for visualization and plotting. 

ANSYS-FLUENT solves the generalized Navier-Stokes equations using the finite volume 

technique. It has a pressure- and a density-based solver to calculate the incompressible and 

compressible flows, respectively. The standard version of ANSYS-FLUENT has a 

comprehensive suite of models to represent conduction, convection, and radiation heat transfer 

with options to simulate phase change and solidification-melting phenomena. Radiation and 

phase change models were used extensively in the simulations described in this report. A 

number of multiphase flow modeling techniques including the volume of flow method, 

Eulerian-Eulerian model, mixture model, and the discrete particle tracking methods are 

available in the solver. The mixture multiphase model was used to account for phase 

dependent quantities in the numerical simulations reported here (ANSYS, Inc., 2009f). 

The standard solver can be customized to meet the requirements unique to a particular 

application by using the user-defined functions (ANSYS, Inc., 2009g). Users can choose from a 

number of turbulence models ranging from zero equation models to large eddy simulation 

techniques available with the standard solver to simulate turbulence as their problems require. 

A number of spatial discretization schemes like the first-order upwind scheme, power law 

scheme, second-order upwind scheme, and central differencing scheme are available with both 

implicit and explicit time integration techniques for temporal terms of the equations. The 

simulation boundary conditions could be defined using the options available with the solver, or 

users can define it through customized functions. Most of the standard boundary conditions 
2-2 

such as the velocity and pressure inlet boundaries, outflows, periodic conditions, and solid wall 

conditions are available with the standard package. A large database of fluid and solid 

properties is also provided with the solver to model the presence of fluids, solids, and mixtures. 

Special boundary conditions and the modifications of equations needed for swirling and rotating 

flows and flows with nonstationary reference frames can be input as user options in the solver. 

ANSYS-FLUENT uses an unstructured grid and supports a number of grid elements such as 

hexahedral, polyhedral, prismatic, and tetrahedral. A combination of hexahedral and tetrahedral 

meshes was used for building the models described in this report. The solver employs the 

Message Passing Interface routines for parallel processing on a number of platforms including 

Microsoft® Windows® NT, UNIX, and different variants of LINUX. The three-dimensional models 

described in Chapter 4 were solved using the parallel processing capability of ANSYS-FLUENT 

and were run on a LINUX cluster at the Center for Nuclear Waste Regulatory Analyses using 

18 compute nodes. 

2.2 Technical Requirements of Customized Functions 

A set of customized functions was developed and plugged into the baseline solver to model the 

water-vapor –liquid phase change process. For the present application, three distinct 

possibilities were identified, where liquid-vapor phase change can take place. 

1. Liquid water available at a relatively hot surface will evaporate from the surface and will 

use the hot surface as the source of latent heat. 

2. Water vapor will condense to form liquid water in contact with a relatively cold surface 

and will use the cold surface as the sink of latent heat. 

3. If the local temperature in certain areas drops below the saturation limit, moist air will 

partially condense to form liquid droplets. Condensation will continue until local 

equilibrium is achieved. In this situation, latent heat of condensation will be released to 

the surrounding air. 

Assumptions were made to analytically quantify the mass and energy exchange as a result of 

these processes and to include them within the baseline simulation (discussed in the following 

section). The analytical formulation was subsequently coded as a User Defined Functions 

(UDFs) and was tagged with the main solver. The effects of these processes were incorporated 

as source terms of the governing equations and suitable boundary conditions. 

2.2.1 Nomenclature 

mሶ HమO ᇱᇱᇱ — Volumetric mass source of water vapor in a computational cell (kg/m3

-s) 

mሶ ୟ୧୰

ᇱᇱᇱ — Volumetric mass source of air in a computational cell (kg/m3

-s) 

mሶ

ᇱᇱᇱ — Total (air + water) volumetric mass source (kg/m3

-s) 

mሶ HమO ᇱᇱ — Mass flux of water vapor through boundary layer (kg/m2

-s) 

mሶ ୟ୧୰

ᇱᇱ — Mass flux of air vapor through boundary layer (kg/m2

-s) 

mሶ

ᇱᇱ — Total (air + water) mass flux through boundary layer (kg/m2

-s) 

n Direction normal to the wall and assumed to be the y direction in the 

 formulation 

ρ — Mixture density (kg/m3

) 
2-3 

D — Diffusion coefficient (m2

/s) 

y — Vertical distance from bottom wall (m) 

v — Normal fluid velocity at the condensed water mixture interface (m/s) 

ωHమO — Mass fraction of water vapor 

ωHమO ሺT୵ୟ୪୪ሻ — Saturation mass fraction of water vapor at wall temperature Twallωୟ୧୰ 

 — mass fraction of air 

mୡ୭୬ୢ — Maximum condensable water vapor mass in volumetric condensation 

m୤,ୡୣ୪୪ — Cell water vapor mass fraction 

m୤,ୱୟ୲୳୰ୟ୲ୣୢ — Cell saturated mass fraction at cell temperature 

f୰ୣ୪ୟ୶ — Relaxation factor for volumetric condensation 

mሶ — Condensation rate due to volumetric condensation 

S — Supersaturation Ratio 

ሺP୴ሻSSL — Vapor Pressure under supersaturation limit (Pa) 

Pୱୟ୲ሺTሻ — Saturation Vapor Pressure at Temperature T (Pa) 

T — Cell Temperature (K) 

J — Net flux of the number of droplet embryos {for droplet nucleation106

 m-3/s 

 [2.831×104

 ft−3/s]} 

σ — Surface tension between water droplets and surrounding moist air (N/m) 

Mഥ — Molecular Weight of moist air (kg/Kmol) 

R — Universal Gas Constant (J/Kmol/K) 

NA — Avogadro Number [6.022141×1023/mol] 

KB — Boltzmann Constant [1.38×10−23 J/K (5.657×10−24 ft-lb/o

R] 

vl — Specific volume of liquid water at cell temperature(m3

/kg) 

Aୡୣ୪୪ ୵ୟ୪୪ — Face area of a computational cell at wall (m2) 

Vୡୣ୪୪ — Volume of a computational cell at wall (m3

) 

2.2.2 Wall Condensation 

The model for wall condensation was adopted from the study of Bell (2003) for condensation of 

superheated water on cold flat plate. It is conceptualized that condensation occurs as water 

vapor diffuses through the mass transfer boundary layer at the surface and on contact with the 

cold surface releases the latent heat to form a liquid film of water. The liquid film attains the wall 

temperature and the vapor diffusion and condensation continues at the liquid film-gas phase 

(air-vapor mixture) interface. This process is shown schematically in Figure 2-1. The rate of 

condensation is determined by the rate of vapor diffusion, which in turn depends on the 

difference in water vapor concentration across the boundary layer. On one end of the boundary 

layer is the bulk mixture and the concentration depends on the bulk flow parameters. The other 

end of the boundary layer is the condensed liquid water film and gas mixture interface. The 

present study assumes a local equilibrium at the liquid film-gaseous air water vapor mixture 

interface, which implies that the water vapor concentration at the interface is equal to the 

saturation value at the film temperature. Therefore, the film interface temperature dictates the 

concentration difference across the boundary layer and in turn controls the condensation rate. 

The following assumptions are made regarding wall condensation in this analysis. 

(i) Condensation occurs due to diffusion of water vapor through the species boundary layer. 

This assumption is used for deriving the mathematical expression for condensed 

water volume. 
2-4 

Figure 2-1. Schematic of the Condensation Model 

(ii) The liquid film is in thermal equilibrium with the boundary wall and has the same 

temperature as the wall. 

(iii) As discussed before, a saturation condition is also assumed at the liquid film-gas (air 

water-vapor mixture) interface. This implies that at the interface the partial pressure of 

water vapor is equal to the saturation vapor pressure at the cold wall temperature. This 

assumption, coupled with assumption (ii), indicates that the wall temperature dictates the 

water vapor concentration at the interface and in turn controls the total vapor diffusion 

and concentration rate across the boundary layer. 

(iv) Any effect of the movement of liquid film is not considered. The liquid film may slide due 

to fluid shear and gravity, but it is assumed that the resulting motion does not affect the 

convecting airflow or subsequent condensation. 

(v) Dropwise condensation is not considered and the analysis is restricted only to 

film condensation. 

(vi) Any thermal resistance offered by the liquid film is not considered. 

(vii) The liquid film consists of water and no other impurities exist. 

(viii) Thermodynamic equilibrium is assumed at the liquid-film-air water-vapor 

mixture interface. 

These assumptions are reasonable for the intended use of the developed UDFs, which is the 

condensation of water vapor within an underground drift. At the same time these assumptions 

also limit the use of the developed UDFs only for specific scenarios. Application of these UDFs 

to simulate a different environment will require a careful review of the effect of these 

assumptions on the problem to be modeled. 
2-5 

As mentioned previously, the condensation rate under these conditions is limited by the 

diffusion of water vapor towards the cold wall. The species mass flux equation for water vapor 

(Bird, et al.,1960; Bell, 2003) and air at the cold wall is written as 

mሶ HమO ᇱᇱ ൌ െρD

∂ωHమO

∂n ൅ ωHమO ρv (2-1)

mሶ ୟ୧୰

ᇱᇱ ൌ െρD

∂ωୟ୧୰

∂n ൅ ωୟ୧୰ ρv (2-2)

From the definition of mass fraction we can conclude that 

ωୟ୧୰ ൌ ሺ1 െ ωHమO ሻ

(2-3)

∂ωHమO

∂n ൌ െ ∂ωୟ୧୰∂n (2-4)

The total mass flow through the boundary layer is given as 

mሶ

ᇱᇱ ൌ mሶ ୟ୧୰

ᇱᇱ ൅ mሶ HమO ᇱᇱ ൌ ρv

(2-5)

The air, however, does not dissolve into the liquid water film, and it is assumed that the mass 

flow rate of air is zero at the interface between liquid water and air-water vapor mixture. This 

assumption, along with Eqs. (2-1), (2-3), and (2-4), leads to the following relationship. 

ρv ൌ 1

൫ωHమO െ 1൯ ρD

∂ωHమO

∂n (2-6)

From the mass continuity of a single computational cell near the boundary, the following 

relationship can be established 

mሶ

ᇱᇱᇱ ൌ ρv

Aୡୣ୪୪ ୵ୟ୪୪

Vୡୣ୪୪ (2-7)

The ρv term on the right hand side of Eq. (2-7) represents the flux that flows into a 

computational cell adjacent to a wall, and the product of the flux with the cell surface area at 

the cell-wall interface provides the total mass flow through the cell. It is assumed that the 

flow variables and properties do not vary considerably within a computational cell adjacent to 

the wall. This assumption is reasonable as long as the wall-adjacent cells are relatively small 

in size. 

The total volumetric source term for the equation can be derived using Eqs. (2-6) and (2-7) as 

mሶ

ᇱᇱᇱ ൌ 1

൫ωHమO െ 1൯ ρD

∂ωHమO

∂n

Aୡୣ୪୪ ୵ୟ୪୪

Vୡୣ୪୪ (2-8)
2-6 

The volumetric source term for water vapor in the cell can be calculated using Eq. (2-9) 

mሶ HమO ᇱᇱᇱ ൌ ωHమO ൈ mሶ

ᇱᇱᇱ (2-9)

The mathematical expressions in Eqs. (2-1) through (2-9) were used in the baseline solver for 

calculating the source terms. 

2.2.2.1 Source Term Specification 

The baseline solver has the option of adding suitable source terms through customized user 

defined functions. The source terms for the governing equations are derived based on 

Eqs. (2-6) and (2-7) using the following method. 

(i) The total volumetric source term obtained from Eq. (2-8) is used in the 

continuity equation. 

(ii) The volumetric source term for water vapor obtained from Eq. (2-7) is used in the 

water-vapor species equation. 

(iii) The momentum source terms for each direction were derived simply by taking the 

product of the total volumetric source term and the velocity in that direction. 

(iv) The energy source term for the air-water vapor mixture was derived by taking the 

product of the total volumetric source term with the water vapor enthalpy value at the 

respective cell temperature. For wall condensation, the energy source term does not 

include any latent heat transfer term, because of the assumption that the liquid film and 

the boundary wall act as the sole energy sink for phase change. However, during 

volumetric condensation, the source term due to phase change is considered. 

(v) For solid bodies present in the system, the energy source term is derived by taking the 

product of the volumetric source term for water vapor and the enthalpy of condensation. 

The total energy gain due to condensation of water vapor is applied to the solid body 

that contains the condensation surface. 

These source terms were calculated through an UDF and were linked with the baseline solver. 

2.2.2.2 Boundary Condition Specification 

The other important aspect of modeling wall condensation is specification of the proper 

boundary condition. As per assumption (iii) in Section 2.2.2, a saturation condition exists at the 

interface between the liquid film and gaseous air-water vapor mixture. Thus, the relative 

humidity can be set to 100 percent and the water vapor pressure set equal to the saturation 

value at surface temperature. The corresponding mass fraction of the water vapor at the wall 

that contains liquid water will be the saturation mass fraction. The boundary condition at the 

wall is specified according to the equation below

ωHమO ൌ ωHమO,ୱୟ୲ ሺT୵ୟ୪୪ሻ

(2-10)
2-7 

Eq. (2-10) accounts for the wall temperature effect on the condensation process. The wall 

temperature dictates the water vapor mass fraction specified at the wall and that in turn affects 

the vapor diffusion process described in Eq. (2-8). The saturation mass fraction of water vapor 

was determined from the ratio of saturation water vapor pressure and total cell pressure. The 

saturation water vapor pressure at cell temperature was determined using the functional 

relationship provided by Popiel and Wojtkowiak (1998). 

The developed UDFs were used to specify the mass fraction profile as a function of the cell 

temperature. This UDF was linked as the water vapor species boundary condition at specific 

walls that have liquid water for evaporation or cold temperature conducive for condensation. 

2.2.3 Wall Evaporation 

The condensation model proposed by Bell (2003) was altered and adapted for evaporation. 

Evaporation took place from surfaces that had a film of water and had a temperature above the 

surroundings. Evaporation took place as the liquid water changed phase and diffused through 

the mass transfer boundary layer from the liquid layer to the air-water vapor mixture. The 

conceptual description of the evaporation process is described in Figure 2-2. Like the 

condensation rate, the evaporation rate is dependent on the rate of water vapor diffusion 

through the boundary layer. The diffusion rate of water vapor is proportional to the difference in 

water vapor concentration across the boundary layer. At the liquid-film-gaseous air water vapor 

mixture interface, an equilibrium condition is assumed, meaning that the concentration of water 

vapor at this interface is equal to the saturation concentration at the film temperature. 

Therefore, the film temperature dictates the rate of water vapor diffusion, which is assumed to 

be the evaporation rate. 

The following assumptions are made regarding wall condensation in this analysis. 

(i) Like condensation, evaporation is a multiphase process involving flow and coupled heat 

and mass transfer in both the liquid film and vapor. The physics can, however, be 

simplified as a scenario dominated by vapor diffusion through the species boundary 

layer. In this study, evaporation modeling assumes that evaporation is limited by the 

diffusion process through the species boundary layer. The mathematical expression for 

evaporated water vapor mass is derived based on this assumption. 

(ii) The liquid film has the same temperature as the boundary wall. 

(iii) The liquid water surface is saturated with vapor. Therefore, the partial pressure of water 

vapor is equal to the saturation vapor pressure at the liquid film temperature. This 

assumption, coupled with assumption (ii), establishes the influence of wall temperature 

on the evaporation rate. 

(iv) A continuous water film is present on the surface where evaporation takes place. 

(v) Motion of this liquid does not affect the evaporation process or the flow of the moist air 

water vapor mixture. 

(vi) Thermal resistance of the liquid film is minimal and is not considered in the analysis. 
2-8 

Figure 2-2. Schematic of the Evaporation Model 

(vii) The liquid film consists of water and no other impurities exist. 

(viii) Thermodynamic equilibrium is assumed at the liquid-film-air water-vapor 

mixture interface. 

As the assumptions indicate, the basic heat and mass transfer processes in evaporation are 

similar to that of condensation. However, evaporation causes a net mass flow of water vapor 

from the surface to the flow through the mass transfer boundary layer. The basic process of 

mass transfer is diffusion controlled and the diffusion is driven by the concentration difference 

between the liquid water film and the main moist airflow. As a result, Eqs. (2-1) through (2-9) 

can be used to determine the volumetric source per computational cell generated due to 

evaporation, but the direction of mass flow needs adjustment. The system loses mass due to 

condensation and gains mass due to evaporation. To reflect this, the mass and water vapor 

species source terms derived in Eq. (2-8) will take the following form for evaporation modeling. 

mሶ

ᇱᇱᇱ ൌ െ1

൫ωHమO െ 1൯ ρD

∂ωHమO

∂n

Aୡୣ୪୪ ୵ୟ୪୪

Vୡୣ୪୪ (2-11)

Equation 2-11 is used to derive source terms for the governing equations for the evaporation 

process. The method of calculating the source terms for governing equations from the 

volumetric mass source term is described in section 2.2.2.1. The boundary condition for 

evaporation was specified using a mass fraction profile was derived based on the saturation 

condition at the wall temperature. As it is assumed that a saturation condition exists at the liquid 

film surface, the relative humidity is fixed at 100 percent. The relative humidity value was used 

to determine the mass fraction of water vapor. Equation (2-12) was used to specify the species 

boundary condition at the wall where evaporation is expected to take place. 

ωHమO ൌ ωHమO,ୱୟ୲ ሺT୵ୟ୪୪ሻ (2-12)
2-9 

The baseline ANSYS-FLUENT solver does not have the option to specify the required species 

boundary condition. Therefore, a UDF was developed based on Eq. (2-12) to implement the 

boundary condition and was linked with the baseline solver. 

2.2.4 Volumetric Condensation 

Volumetric condensation takes place if convecting moist air encounters a local temperature that 

is below the saturation temperature. As the fluid temperature drops, relative humidity increases 

gradually and reaches 100 percent at the saturation temperature. A further drop in temperature 

will result in condensate formation. A number of assumptions are made with regard to 

volumetric condensation modeling. 

(i) Condensed liquid water forms a mist and is convected with the moist air as a 

separate phase. 

(ii) Volumetric condensation formation is assumed to take place only due to droplet 

nucleation; bubble nucleation is not considered in the analysis. 

(iii) Coalescence, breakup or subsequent evaporation of the liquid droplets is 

not considered. 

(iv) Wall interaction with the liquid droplets is not taken into account. This implies the liquid 

droplets behave the same way as moist air when it comes in contact with a wall. 

(v) A mixture model is used to solve the governing equations for the two-phase ‘moist’ 

air-liquid droplet flows. Specifically, the model solves the momentum, continuity, and 

energy equations for the mixture, the volume fraction equations for the secondary 

phases, and algebraic expressions for the relative velocities. 

Based on these assumptions, two different approaches were adopted for modeling the 

volumetric condensation process and are described in the next section. These two approaches 

are called equilibrium and nonequilibrium models, based on their treatment of the 

thermodynamic condition in a computational cell. The equilibrium approach was used by 

Hijikata and Mori (1973) and Brouwers and Chester (1992a,b) whose formulation depended on 

boundary layer analysis. Later Fox, et al. (1997) adopted this approach to study condensation 

and mist formation in closed enclosures. Fox, et al. (1997) also formulated a nonequilibrium 

model and a modified critical saturation model and showed that the experimental data lie 

somewhere between the equilibrium and nonequilibrium predictions. Kang and Kim (1999) 

used both equilibrium and nonequilibrium methods to characterize film condensation of a 

supersaturated steam-air mixture on a flat plate. They found the model performance was 

dependent on the temperature difference between the supersaturated hot gas and cold 

condensation plate. At a low temperature difference, the computed results matched 

experimental observations, with no mist generation. On the other hand, at a high temperature 

difference, the nonequilibrium model performed better. 

2.2.4.1 Equilibrium Model 

In this approach, it is assumed that the condensation process will continue until local 

thermodynamic equilibrium is reached. This implies that supersaturation of the mixture is not 
2-10 

allowed and any water vapor that is in excess of the saturation limit will condense. Based on 

this understanding, the total mass available for condensation in a single computational cell can 

be expressed as 

mୡ୭୬ୢ ൌ ρሺm୤,ୡୣ୪୪ െ m୤,ୱୟ୲୳୰ୟ୲ୣୢሻ (2-13)

Where mf,saturated can be determined based on the ratio of cell saturation vapor pressure and cell 

total pressure. The cell saturation vapor pressure was determined using the functional 

relationship provided by Popiel and Wojtkowiak (1998). All of the simulations described for the 

current report assume a steady state behavior, and all the water vapor available for 

condensation will change phase to form liquid as the solver converges to a steady solution. 

Numerical stability considerations, however, limit the quantity of mass that can be condensed in 

a single iteration as the solution progresses. As a result, a relaxation factor was introduced to 

achieve numerical stability. The modified equation provides a condensation rate estimate 

mሶ ൌ f୰ୣ୪ୟ୶ ൈ ρሺm୤,ୡୣ୪୪ െ m୤,ୱୟ୲୳୰ୟ୲ୣୢሻ (2-14)

The relaxation factor f୰ୣ୪ୟ୶ was determined using numerical experiments as it depends on the 

local thermodynamic condition, moisture content and other numerical and modeling parameters. 

It was found that the value of this factor can range between 0.3–0.9. The factor is allowed to 

have high values for a system with lower level of supersaturation. 

2.2.4.2 Nonequilibrium Model 

The nonequilibium model assumes that a level of supersaturation will exist in fluid. In this case, 

not all of the available water vapor mass that is in excess of the saturation limit will condense. 

The limit of supersaturation is determined using kinetic theory and is given by the following 

equation (Carey, 2007) 

S ൌ ሺP୴ሻSSL

Pୱୟ୲ሺTሻ ൌ exp ቈ Eכ ሺെ ln Jכଵ/ଶሻ

2ሺEכሻଵ/ଶ ൅ ሺെ ln Jכଷ/ଶሻ

቉

ଵ/ଶ (2-15)

Expressions of E* and J* are provided in the following equations 

J

כ ൌ MഥJ

NAv୪

ቆ ttMഥ

2σNA

ቇ

ଵ/ଶ

൤ RT

Pୱୟ୲ሺTሻ൨

(2-16)

Eכ ൌ 16πσଷv୪

ଶ

3kBRଶTଷ

(2-17)

The condensation rate calculated using the nonequilibium model is given by the 

following relation. 

mሶ ൌ f୰ୣ୪ୟ୶ ൈ ρሺm୤,ୡୣ୪୪ െSൈm୤,ୱୟ୲୳୰ୟ୲ୣୢሻ (2-18)

Hence the vapor mass available for condensation is reduced by the supersaturation limit S. 
2-11 

2.2.4.3 Source-Term Calculation 

A UDF was developed to calculate the mass transfer rate from the vapor to liquid phase as 

described in Eqs. (2-14) and (2-18) (Carey, 2007; Fox et al. 1997). Based on these quantities, 

the source terms were calculated for species, mass, momentum and energy equations. For the 

species and mass conservation equation, the source is obtained directly from Eqs. (2-14) and 

(2-18). The source terms for the momentum equations are obtained by taking the product of 

mass source term with velocity. The energy equation source term was determined by taking the 

product of the mass source and latent heat of condensation. The latent heat of condensation 

for the system pressure was obtained from a functional relationship provided by Popiel and 

Wojtkowiak (1998). 

2.3 General Numerics 

This section describes the numerical parameters that were used in all the simulations 

described in this report. A variety of spatial and temporal discretization schemes as well as 

turbulence models are also available in ANSYS-FLUENT. Simulation parameters were 

chosen as suitable for natural convection flows. For the present simulations, the 

pressure-based solver was chosen, as the flow can be considered to be in the incompressible 

fluid. The two-dimensional Navier Stokes equations were solved using an implicit approach. 

The Semi-Implicit Pressure Linked Equations–Consistent (SIMPLEC) algorithm (Van Doormal 

and Raithby, 1984) was used to treat pressure velocity coupling for numerical stability. The 

third order Monotone Upstream-Centered Schemes for Conservation Laws (MUSCL) were used 

to derive the surface values of different variables for the spatial discretization, used to compute 

the convective fluxes. The Rhie-Chow scheme (Rhie and Chow, 1983) was selected to 

interpolate pressure at the control volume face and to satisfy mass conservation without 

pressure oscillations. Wherever applicable, the shear stress transport k-ω model (Menter, 

1994) was used to simulate turbulence. 

A number of radiation models are available in the baseline CFD solver, including the 

surface-to-surface, Rosseland, Discrete Transfer of Ray and this discrete ordinate model. 

Based on previous studies (Das, et al, 2007, 2008) and guidelines provided by ANSYS-Fluent 

(ANSYS, Inc., 2009a) for modeling radiation in a nonparticipating media with natural convection, 

the discrete ordinate model was chosen for the present study. In this approach, the radiative 

transfer equation for an absorbing, emitting, and scattering medium is solved for a finite number 

of discrete solid angles. In the current study, four angular discretizations were used in each 

direction of the spherical coordinates system. 

2.4 Verification Calculation 

A simulation was carried out to verify the accuracy of the developed UDFs. The 

simulation was adopted from the study of Bell (2003). The study of Bell (2003) simulated 

flow of supersaturated water vapor over a cold flat plate and wall film condensation. 

The numerical study presented in the current report considers volumetric condensation in 

addition to wall film condensation. 

The computational domain, grid, and boundary conditions for the condensation problem are 

shown Figure 2-3. Moist air flows inside the domain at inlet. The lower wall is maintained at a 
2-12 

Figure 2-3. Schematic for Condensation of Humid Air Over Flat Plate 

lower temperature compared to the moist airflow. As a result, condensation takes place at the 

cold wall. Sparrow, et al. (1967) derived an analytical solution for this problem. 

2.4.1 Description of the Numerical Model 

The test case is modeled as a two-dimensional steady laminar two-phase flow problem, and the 

incompressible Navier-Stokes equations are solved without using any turbulence models. 

Species transport equations with inlet diffusion and diffusion energy source are solved for the 

vapor phase. Two different velocities with a fixed water vapor concentration at the inlet are 

studied, and results are compared with experimental data. The fluid in the domain is specified 

as a mixture of water vapor and air, and the density of the fluid is determined using the 

volume-weighted mixing law. The customized UDF to incorporate the source terms is compiled 

and linked with the mass, momentum, energy, and species transport equations using the 

standard ANSYS-FLUENT interface. 

A mixture of water vapor and air is specified as the working fluid in the domain. The free stream 

velocity is specified at the velocity inlet of the domain along with temperature and mass fraction 

of water vapor. A pressure outlet boundary condition with specified backflow temperature and 

water vapor mass fraction is specified at the downstream outlet as well as at the top of the 

domain. A number of customized source terms for the continuity, momentum, energy, and 

species equations are introduced through UDFs to model the effect of mass removal due to 

condensation at the cold bottom wall. To capture the boundary layer and the condensation 

process, grids are clustered near the bottom wall. 

Both equilibrium and nonequilibrium models were used separately for the problem. Figures 2-4 

through 2-7 show results obtained using the nonequilibrium model, whereas Figures 2-8 through 

2-11 highlight results obtained from the equilibrium model that accounts for interphase mass 

transfer due to volumetric condensation. Both these test cases were run for two different 

mixture inlet velocities of 1 m/s [3.28 ft/s] and 0.1 m/s [0.33 ft/s] but the mass fraction of vapor in 

the mixture was fixed at 0.47967. 
2-13 

Figure 2-4. Air Mass Fraction Variation in the Flow Field 

Figure 2-5. Comparison of Analytical and Computed Results for Condensation Mass Flux 

for Inlet Velocity = 0.1 m/s [0.33 ft/s] Using Nonequilibium Model 

[1 m=3.28 ft; 1 kg/m2

-s =0.205 lb/ft2-s] 

Figure 2-6. Comparison of Analytical and Computed Results for Condensation Mass Flux 

for Inlet Velocity = 1 m/s [3.28 ft/s] Using Nonequilibium Model 

[1 m=3.28 ft; 1 kg/m2

-s =0.205 lb/ft2-s] 
2-14 

Figure 2-7. Relative Humidity Contours for Simulation Using 

Nonequilibium Model and Inlet Velocity = 1 m/s [3.28 ft/s]

Figure 2-8. Air Mass Fraction Variation in the Flow Field With Multiphase Flow 

and Inlet Velocity = 1m/s [3.28 ft/s]

Figure 2-9. Comparison of Analytical and Computed Results for Condensation Mass Flux 

for Inlet Velocity = 0.1 m/s [0.33 ft/s] Using Equilibrium Model 

[1 m=3.28 ft; 1 kg/m2

-s =0.205 lb/ft2-s] 
2-15 

Figure 2-10. Comparison of Analytical and Computed Results for Condensation 

Mass Flux for Inlet Velocity = 1 m/s [3.28 ft/s] Using Equilibrium Model 

[1 m=3.28 ft; 1 kg/m2

-s =0.205 lb/ft2-s] 

Figure 2-11. Relative Humidity Contours for Simulation With Multiphase Species 

Transport and Inlet Velocity = 1 m/s [3.28 ft/s] 

2.4.2 Test Results for Nonequilibrium Model 

The air mass fraction contours for an inlet velocity of 1 m/s [3.28 ft/s] are shown in Figure 2-4. 

The thickness of the air mass fraction increases along the cold wall because the mixture loses 

water vapor due to condensation and the mass fraction of water vapor in the mixture decreases. 

Consequently, the mass fraction of air in the mixture increases near the bottom wall of the 

domain. This is consistent with the understanding of the physics of film condensation on cold 

flat plates. 

Figure 2-5 compares condensed mass flux at the cold bottom wall for an inlet velocity of 0.1 m/s 

[0.33 ft/s]. In the downstream region, the pattern of the computed results is in good agreement 

with the analytical solution, though the computed result slightly under predicts the data. 
2-16 

The analytical and computed solutions have some deviation near the leading edge of the plate 

{0–0.1 m [0–0.328 ft]}, though near the trailing edge {0.8-1 m [2.62–3.25 ft]} the deviations do 

not exceed 10 percent of the overall range. The computed solutions do not take into account 

the boundary layer development and assume a fully developed boundary layer from the leading 

edge. Thus, the leading edge results should be excluded from the comparison and the study is 

considered validated only for the downstream flow. The deviation in the trailing edge occurs 

because the analytical solution makes certain assumptions regarding the boundary layer 

thickness near the plate leading edge that are different from the simulated case. 

Figure 2-6 compares the computed and analytical solution for an inlet velocity of 1 m/s [3.28 ft/s] 

and shows the same trend as in the previous test case with a different inlet velocity. The results 

show some deviation between simulated and analytical solution is within 10 percent based on 

the total range. 

Figure 2-7 shows the contours of relative humidity of the entire flow domain for an inlet velocity 

of 1 m/s [3.28 ft/s]. The results show that though the majority of the domain is either at a 

saturated or unsaturated condition, there is a thick layer near the condensation zone that shows 

supersaturation. As mentioned in Section 2.2.4.2, the nonequilibrium solution does not consider 

complete volumetric condensation where the excess water vapor that causes supersaturation is 

partially allowed to remain in vapor phase. As a result, a small section of fluid near the 

boundary layer is supersaturated. 

2.4.3 Test Results for Equilibrium Model 

Figure 2-8 shows calculated concentration of air in the domain using equilibrium model. The 

result is similar to the concentration distribution obtained using nonequilibrium model. This 

similarity is due to the lack of volumetric coordination in the main flow away from the wall. As 

there is negligible volumetric condensation, equilibirum and nonequilibrium models work the 

same way and calculates similar concentration distribution. Due to a high wall condensation 

rate, the air mass fraction near the bottom cold wall is affected. Inside the domain, away from 

the wall, volumetric condensation is the only mechanism that can affect species distribution. 

Because the volumetric condensation rate is orders of magnitude less than the wall 

condensation rate, the species distribution appears to be unaffected. 

Figures 2-9 and 2-10 show the wall condensation mass flux for two different equilibrium model 

test runs. In general, both the runs reasonably predict the experimental data. 

However, unlike the test cases described in Section 2.4.2 in connection with the nonequilibrium 

model where the computed results slightly under-predicted the condensation rate, here the 

computed solution slightly over-predicts it. The disparity between the computed and the 

analytical solution in the upstream region is higher compared to the nonequilibrium model 

solution. This can be attributed to a number of factors. 

1. The analytical solution was derived for wall condensation only and did not consider 

volumetric condensation. Thus, some disagreement between the equilibrium 

model solution, which considers volumetric condensation, with analytical solution 

was expected. 
2-17 

2. Removal of water vapor from the mixture phase may have caused a steeper gradient of 

species concentration, causing a higher diffusion flux of water vapor towards the wall. 

3. The volumetric condensation process affects the species distribution of water vapor near 

the cold wall and can affect the diffusion flux of water vapor through the boundary layer 

that ultimately affects the condensation rate. 

4. A number of assumptions regarding boundary layer development were made in the 

analytical solution near the flat plate leading edge which are not present in the 

computational solution. The boundary layer development in the computational solution 

did not follow the prescribed profile specified for the analytical solution. 

Figure 2-11 shows the relative humidity contours for the flow domain. The results show that 

the relative humidity pattern obtained using the equilibrium model are different from that 

obtained using the nonequilibium model. In the results obtained using equilibrium model, 

the supersaturation layer has been reduced significantly. Ideally, the entire domain should 

either be at an unsaturated or saturated condition. The nonequilibium model was not able 

to completely remove the extra water vapor responsible for saturation due to stability issues 

with the numerical solution. This problem is encountered for systems with a higher level of 

moisture content and will not pose a problem for systems with moderate or low levels of 

moisture content. 

The verification calculation in this chapter shows that the UDFs for volumetric and surface 

condensation were able to calculate a flow field that matched experimental data with less than 

10 percent deviation. The nonequilibrium model produced relatively better results than the 

analytical solution which did not consider volumetric condensation. The relative humidity 

calculated using the equilibrium model was, however, more realistic as supersaturation 

was not allowed. 

2.5 Summary 

This chapter discusses the numerical approach utilized to study the condensation-evaporation 

cycle along with natural condensation in the cold-trap process. ANSYS-FLUENT Version 12.1 

is used as the general baseline solver for the problem. A set of customized UDFs were 

developed to model the condensation and evaporation processes. The mathematical bases for 

the wall evaporation and condensation models were discussed. The volumetric evaporation 

process was modeled using two different approaches depending on the treatment of 

thermodynamic equilibrium of the system. The equilibrium model assumed that the volumetric 

condensation process will continue until a saturation condition is attained. On the other hand, 

the nonequilibrium model allowed partial supersaturation and the level of supersaturation was 

obtained using kinetic theory. A verification calculation was done where an analytical solution 

for the condensation rate of hot moist air on a cold plate was compared with simulated data. 

This study provides confidence in the analytical formulation of the condensation and 

evaporation processes and the implementation of these mathematical expressions as UDFs. 
3-1 

3 CONDENSATION CELL EXPERIMENT AND VALIDATION 

Two detailed laboratory tests were conducted to verify the overall modeling capabilities of the 

ANSYS-FLUENT® solver in predicting the flow and thermal field. The first set of experiments 

used a small-scale laboratory model to validate the user defined modules for calculating the 

wall and volumetric evaporation-condensation process. The environment within the test 

chamber was well controlled and regulated. The experiments were conducted to obtain 

measurements of water transport rates in a simple two-dimensional natural convection flow. To 

simplify the validation experiments, low temperatures and temperature gradients were used so 

that radiation could be neglected. The intent of these experiments was to make a quantitative 

comparison between the experimentally obtained and computed condensation rates and air 

temperatures. This chapter will first describe the small-scale experiments and their results for 

temperature distributions and condensation rates. Subsequently, development of the 

corresponding numerical model will be described and results from the computational fluid 

dynamics (CFD) analysis will be compared to the measured water transport rates. 

3.1 Description of Water Transport Experiments 

Laboratory tests were conducted to measure the amount of water transported from a heated 

water source to a cold (condensation) plate. Multiple tests were conducted by varying the heat 

rate to the water source and varying the temperature of the cold plate. Once steady state 

temperatures were attained, water transport rates were determined by collecting condensation 

draining from the cold plate. A description of the test chamber is given below. 

3.1.1 Test Chamber 

A schematic diagram and photographs of the test chamber are shown in Figures 3-1 and 3-2, 

respectively. The enclosure was made of polycarbonate sheet, acrylic sheet, and aluminum. 

Water evaporates from the heated tray shown on the bottom left of the enclosure. Water was 

continuously added to the water tray to maintain the level. Water was condensed and collected 

on the cooled aluminum plate shown on the right side of the test enclosure. Instrumentation 

was included to measure the evaporator (water) temperature, condenser temperature, air 

temperatures, and the water condensation rate. 

Figure 3-1. Side View Schematic of Laboratory Test Chamber 

Warm Water

Water

Chiller

Cool

Plate

Electric Heater

T

T Water

Supply

Condensation

Collection
3-2 

Figure 3-2. Side View Photographaboratory Test Chamber 

The chamber was 0.584 m [23 in] long, 0.154 m [6 in] tall and 0.305 m [12 in] deep. The tray 

holding the source water was 0.076 m [3 in] long and ran the entire 0.584 m [12 in] depth of the 

chamber. The top wall of the chamber was made from 0.025 m [1 in] thick acrylic. The left end 

wall and bottom wall were made of 0.025 m [1 in] thick polycarbonate. The water tray placed at 

one end of the bottom wall, as shown in Figure 3-2, was made of aluminum. The front and back 

walls of the chamber were made of 0.013 m [0.5 in] thick polycarbonate. The cold plate was 

made of 0.025 m [1 in] thick aluminum and had channels drilled through it to allow the 

circulation of cooling fluid through the plate. Though not shown in Figure 3-2, the entire 

chamber was covered with a 0.025 m [1 in], and later 0.050 m [2 in], thick shell of polystyrene 

foam insulation to reduce heat transfer between the experimental chamber and the 

surroundings. The chamber was placed on top of a table and to reduce heat transfer, two 

sheets of polystyrene foam insulation were placed between the chamber and the table top. 

External thermal energy was provided for warming the water using an electric heater, which was 

attached directly to the bottom of the aluminum water tray. The heater pad had the same 

dimensions as the water tray, which was 0.076 m [3 in] long by 0.305 m [12 in] deep. The 

power to the heater pad was supplied by a variable power source. The goal of the heating and 

cooling component design in the chamber was to create a two-dimensional convection cell in 

the plane along the length and height of the chamber, with minimal variation in the third 

direction. Temperature measurements were made at the centerline of the 0.305 m [12-in] depth 

of the rectangular domain. 

Thermocouples were installed in the test chamber to measure the temperature and their 

locations are shown in Figure 3-3. The thermocouples were Omega, Type K, model number 

5TC TT K 30 72 with a bare wire junction and 25 ×10−5

 m [0.01 in] diameter wire. A calibration 

check showed that all of the thermocouples were within 0.2 °

C [0.36 °F] of the standard 

thermometer reading. 

Water Tray Cold Wall 
3-3 

Figure 3-3. Top and Side Views Schematic of Laboratory Test Chamber. 

All Dimensions Are in Inches. T1 to T9 Are Thermocouples. [1 m = 39.37 in] 

Water that condensed on the cold plate was drained off and collected in a graduated cylinder to 

measure the condensation rate. Depending upon the condensation rate and collection time, 

different size graduated cylinders were used to collect the fluid to improve measurement 

resolution. Several tests were initially performed without water to test the estimation of 

convection velocities and temperatures in a simpler regime and prepare the system for actual 

tests. Results of these initial runs are not discussed in this report. 

3.1.2 Observations and Test Measurements 

Several tests were conducted to measure the moisture transport rate from the hot water tray at 

the bottom wall at one end of the chamber to the cold plate at the sidewall of the chamber. For 

each thermal loading and cold plate temperature, the water condensation rate at the cold plate 

was recorded after steady state was established within the experimental enclosure. The test 

results are shown in Table 3-1. The values given in Table 3-1 are averages measured over 

several hours after steady state temperatures were reached. 

Light condensation was observed to form on the inside surface of the chamber along the top 

wall and the sidewalls after the initial set of tests (Tests 1–3). Starting from the edge of the 

cold-wall, the top-wall surface covered by condensate spanned a distance of about 0.305 m 

[12 in]. Closer to the cold plate, higher condensate mass and larger condensate droplet size 

Side

View

Top

View

Water Trough

4.0

11.5

4.0

8.0

0.1

6.0

2.0

T1

T2

T3,T4,T5, T6,T7 T8,T9

1.5

5.0 1.0 4.5

3.0

Acrylic Sheet, 1" Thick

Styrofoam Sheet, 1" Thick

T1,T2

T9

T8

T5

T4

T3

T7

T6

Cold

Plate

End

Wall

Water Tray

Walls 
3-4 

Table 3-1. Test Results and Test Conditions. 

Temperature (Temp.) Sensor Locations (T3–T9) Are Shown in Figure 3-3. 

Test 

No. 

Condensation 

Rate (ml/hr) 

Water 

Temp. 

(C ) 

Cold 

Plate 

Temp. 

(C ) 

Air 

TempLower 

(C) 

Air 

TempMiddle 

(C) 

Air 

TempUpper 

(C) 

Temp. Top 

Chamber 

(C) 

Temp. 

Top 

Insulation 

(C ) 

T1, T2 T6, T7 T3 T4 T5 T8 T9 

1 8.8 38.6 10.7 22.9 26.0 27.2 25.4 23.7 

2 16.8 46.1 10.8 26.6 30.2 31.5 29.0 24.1 

3 25.6 54.3 10.9 30.9 35.2 36.3 33.1 24.5 

4 7.2 32.3 5.3 18.8 21.6 23.2 22.4 24.5 

5 12.1 39.4 5.4 21.5 24.9 26.6 24.7 23.1 

6 20.9 47.7 5.5 25.1 29.1 30.9 29.0 23.7 

7 4.7 25.3 5.2 15.6 17.9 19.2 18.7 23.3 

8 2.5 19.5 5.2 14.3 16.3 17.4 18.2 24.2 

9 1.7 26.2 19.1 22.6 23.5 24.1 24.0 25.4 

10 6.6 38.2 19.1 26.3 28.3 29.5 27.9 25.1 

11 20.8 50.9 19.3 32.1 35.3 36.8 34.4 25.1 

12 12.8 45.0 19.2 29.2 32.1 33.4 31.5 26.2 

13 8.1 46.1 29.5 33.7 35.5 36.6 34.8 25.1 

14 4.2 39.7 29.4 31.3 32.2 33.1 31.5 23.9 

15 0.5 34.0 29.4 29.7 30.1 30.5 29.6 25.2 

was observed. Condensate mass accumulation was also observed on the side walls but the 

condensation layer did not reach all the way to the bottom of the chamber. A rectangular 

strip, approximately 0.051 m [2 in] high along the interface between the bottom and side wall, 

was free of condensation. As condensate formation was observed on the side walls of the 

chamber in Tests 1–3, an analysis was performed to determine the effect the heat loss from 

the chamber on condensation. An analysis of the energy balance was performed considering 

the heater power and the heat loss to the ambient environment. The heat loss through the 

insulation was estimated using the insulation thickness, chamber surface area, insulation 

thermal conductivity 0.026 W/m K [0.15 BTU/hr ft °F], and the temperature difference measured 

across the insulation on the top of the chamber. The amount of water that could condense in 

the chamber from the heat lost to the ambient environment was estimated from the heat loss 

rate and the heat of vaporization for water, estimated to be 0.63 W hr/g at 100 °C [975 BTU/lb 

at 212 °F]. The heat loss rate for Test 3 was estimated to be 5.4 W [18.426 BTU/hr]. This 

heat loss rate corresponds to a condensation rate of 8.6 ml/hr [0.002 gal/hr] of water. For this 

test, the measured condensation rate on the cold wall was 25.6 ml/hr [0.007 gal/hr]. Hence the 

side wall condensation was approximately 29 percent of the cold wall condensation and cannot 

be neglected. The simplified heat loss analysis, however, neglects any cooling of the water or 

air in the chamber and assumes all heat loss to the ambient environment results in 

condensation. The analysis revealed that a significant amount of heat loss was occurring and to 

reduce it, after Test 3, the shell of insulation surrounding the test chamber was increased to a 

0.051-m [2-in] thickness. 

Data from Tests 4-15 in Table 3-1 are summarized in Figure 3-4, where the variation of the 

condensation rate is plotted as a function of the temperature difference between the hot water 

source and the cooled condensation wall. Four curves, labeled as 5, 10, 20, and 30 °C 

[41, 50, 68, and 86 °F] curves, represent different absolute temperature regimes in the cold wall 

plate. Figure 3-4 shows that condensation increases both with the temperature difference 
3-5 

Figure 3-4. Experimental Results Showing the Change in the Condensation Rate As a 

Function of the Difference Between the Heated Water and Cold Plate Temperatures 

 

between the hot source water and cold condensation plate and also with the absolute temperate 

regime of the experiment. 

3.2 Numerical Modeling of the Moisture Transport Experiments 

A numerical model was developed using the commercial CFD package ANSYS-FLUENT 

Version 12.1 to replicate the experimental study. This general purpose package uses a 

control-volume-based technique to convert a scalar transport equation to an algebraic equation 

that is solved numerically. A general description of the package is provided in Chapter 2. 

Simulation parameters were specifically chosen to suit natural convection flows. The 

experimental setup was designed to have a single condensation-evaporation cell in a 

two-dimensional plane with almost no variation in the third dimension. To ensure the 

two-dimensionality of the flow field, no heat source or sink was placed on the side walls and 

as mentioned in the previous section, heat loss from the side walls was minimized. Hence a 

two-dimensional domain was used in the computational study. Rayleigh number calculation 

indicates that the flow could be transitional or turbulent, depending on the temperature 

difference between the hot water and cold plate. Hence, resolving turbulence quantities in the 

flow was deemed important and the two-equation shear stress transport k-ω turbulence model 

was used in the study. 

3.2.1 Numerical Method 

General mass transfer modeling techniques used in this study were described in Chapter 2. 

The experimental study revealed that there were three distinct mass transfer and transport 

processes that took place within the chamber. 

(i) Mass transfer at the walls due to condensation and evaporation that is modeled using 

suitable boundary conditions 
3-6 

(ii) Transport of evaporated water vapor from the hot source location to the cold wall due to 

buoyancy driven natural convection 

(iii) Condensation of water vapor within the flow domain, where the temperature locally 

drops below the saturation point 

The third process would generate liquid water that will then be transported with the air and water 

vapor mixture. The quantity of condensed liquid water would be very small, and the effect of 

this liquid phase on overall gas velocity will be relatively small. The thermodynamic effect, 

however, could be significant due to latent heat exchange between phases. As mentioned in 

Chapter 2, the mass transfer at the wall due to evaporation and condensation is modeled using 

suitable boundary conditions that assume an equilibrium condition at the liquid-water-gaseous 

mixture interface. This approach also assumes only film condensation occurs (i.e., droplet 

condensation is excluded). Also, any phase change heat transfer takes place within the wall 

and does not affect the gas phase. The condensation or evaporation rate is calculated using 

the water vapor diffusion rate. The volumetric condensation in the domain was simulated using 

both equilibrium and nonequilibium techniques. The theoretical basis of these techniques and 

their formulation were described in Chapter 2. The mass transfer boundary conditions at the 

walls were coded as User Defined Functions and loaded into the solver as surface profiles. The 

source terms due to mass exchange between phases were added to the mass, momentum and 

energy equations. 

3.2.2 Modeling Results 

Figure 3-5 shows the computational grid domain and the boundary conditions of the model. A 

total of 17,000 nodes were used in the calculation and this grid dimension was determined using 

a separate grid independence study. It was found that approximately 1.5 × 104

 nodes in the 

domain provided a stable flow field that did not change with grid dimension. Grids were 

clustered near the wall to capture flow, temperature, and concentration gradients. As discussed 

in the previous section, a saturation condition was specified at the water source and the cold 

wall corresponding to the experimental temperature at those locations. The rest of the wall was 

assigned an adiabatic boundary condition as the heat transfer through those walls was 

assumed to be negligible. 

Simulations were performed for the conditions corresponding to Tests 1–14. The general flow 

and temperature field data presented here for both the models correspond to Test No. 11 of 

Table 3-I. 

Figure 3-5. Computational Grid and Boundary Condition 
3-7 

Figures 3-6 and 3-7 show the velocity magnitude contours obtained using the nonequilibrium 

and equilibrium models, respectively. In general, a strong plume is generated from the hot 

water source that moves upward due to the density difference between the hot and the cold 

fluid. This plume reaches the upper wall and then moves toward the cold wall. Another high 

velocity region can be noticed near the cold wall, where the fluid undergoes mass transfer and 

the colder fluid flows back to the water source creating a cyclic convection pattern that encircles 

the whole domain. Though there is similarity in the velocity distribution pattern predicted by the 

two models, the velocity magnitude within the circulation loop calculated by the equilibrium 

model is higher compared to the nonequilibrium model. This may be attributed to the fact that 

volumetric condensation has caused a volume deficit in the enclosure that results in a higher 

velocity of the gas phase. 

Figures 3-8 and 3-9 show the temperature field calculated using the equilibrium and 

nonequilibrium models. A high temperature zone can be identified near the hot water source 

and the temperature adjacent to the cold wall is lower than the rest of the domain. Though the 

temperature patterns predicted by both models are qualitatively similar, some difference in 

magnitude can be noted near the bottom wall of the enclosure. Any difference in the 

Figure 3-6. Velocity (m/s) Contours Using Nonequilibrium Model 

(Without Volumetric Condensation) [1 m/s = 3.28 ft/s]

Figure 3-7. Velocity (m/s) Contours Using Equilibrium Model 

(With Volumetric Condensation) [1 m/s = 3.28 ft/s]
3-8 

Figure 3-8. Temperature (K) Contours Using Nonequilibrium Model 

(Without Volumetric Condensation) [°F = 1.8 °K × −459.4]

Figure 3-9. Temperature (K) Contours Using Equilibrium Model 

(With Volumetric Condensation) [°F = 1.8 °K × −459.4]

temperature distribution pattern is solely due to volumetric condensation. The wall 

evaporation-condensation does not affect fluid temperature inside the domain as latent 

heat exchange at the wall takes place with the solid material of the enclosure. In general, 

the temperature distribution pattern is not greatly affected by volumetric condensate formation 

as the mass of condensate is relatively low. Volumetric condensation, however, can affect 

localized temperature distribution. 

Figures 3-10 and 3-11 illustrate the relative humidity contours of the nonequilibrium and 

equilibrium models. Note that a substantial part of the domain has relative humidity 

slightly above 100 percent for simulations that use the nonequilibrium model. This is 

because supersaturation was allowed in the domain and any possible mass transfer from 

the air-water-vapor mixture was not considered in the calculation. The equilibrium model 

shows a different pattern of relative humidity distribution, where the maximum value is capped 

at 100 percent as any water vapor that contributes to supersaturation condenses to form 

liquid water. Though supersaturation does not normally occur in an enclosed space, 

researchers (Brouwers and Chester, 1992a,b; Kang and Kim, 1999; Carey, 2007) have 

indicated that it can be attained under limited circumstances and the maximum level of 

supersaturation can be determined using kinetic theory. No experimental measurement was 

taken for relative humidity in the chamber. As a result, no conclusion on the accuracy of the 

computed relative humidity pattern using either of the approaches can be made. 
3-9 

Figure 3-10. Relative Humidity Contours Using Nonequilibrium Model 

(Without Volumetric Condensation)

Figure 3-11. Relative Humidity Contours Using Equilibrium Model 

(With Volumetric Condensation)

The nonequilibium and equilibrium models, however, provide two bounding scenarios as these 

models consider the minimum and maximum possible water vapor mass that could undergo 

volumetric condensation. 

Figure 3-12 shows the interphase volumetric mass transfer rate for the equilibrium model. The 

overall volumeric mass transfer rate is small and values are largest in a small region near the 

hot water source and the cold wall. These regions correspond to high relative humidity zones 

obtained using the equilibrium model as shown in Figure 3-11. The mass transfer is expected 

to be maximum in regions where the degree of supersaturation is high. This is because the 

equilibrium model assumes that mass transfer rate is proportional to the difference between the 

supersaturation and saturation water vapor concentrations. 

3.2.3 Comparison Between and Experimental and Computed Data 

Figures 3-13 and 3-14 show the comparison between experimentally obtained and computed 

condensation rates using the nonequilibrium and equilibrium models, respectively. The data are 

grouped together in four distinct sets that correspond to different convection regimes mentioned 

in Section 3.1.2. Generally, both computed data sets matched the trend of the experimental 

data, where the condensate rate increases with increased temperature difference. 
3-10 

Figure 3-12. Interphase Mass Transfer (kg/m2

-s) [1 kg/m2-s = 0.205 lb/ft2-s] 

Due to Volumetric Condensation Using Multiphase Flow

Figure 3-13. Comparison of Experimental and Computed Values of Condensation Rate 

Using Nonequilibrium Model (Without Volumetric Condensation) 

[1gal/hr = 3785.41 ml/hr; °F = 1.8 × °C + 32] 

Though both models overpredict condensation rate, a closer examination of the results shows 

that the computed results using the equilibrium model match the experimental data better, 

especially at a lower temperature difference. For example, the lower water vapor content in the 

gaseous mixture near the cold wall (using the equilibrium model) is lower than the 

corresponding values obtained using nonequilibium model. This is because some of the water 

vapor available within the mixture has already been allowed to condense and form liquid mist in 

the equilibrium model. On the other hand, the entire water vapor content is convected from the 

hot water pool to the cold wall in the nonequilibium model calculations resulting in a higher wall 

condensation rate. 

Figures 3-13 and 3-14 also establish the effect of the temperature difference between the cold 

and hot wall in total mass transfer rate across the domain. The quantity of water evaporating at 

the hot surface depends on the temperature difference between the hot water pan and the bulk 
3-11 

Figure 3-14. Comparison of Experimental and Computed Values of Condensation Rate 

Using Equilibrium Model (With Volumetric Condensation) 

[1gal/hr = 3785.41 ml/hr; °F = 1.8 × °C + 32]

temperature. On the other hand, the condensation rate depends on the temperature difference 

between the cold plate and the bulk fluid (air). Under steady state conditions, the rate of 

condensation will be equal to the rate of evaporation. In essence, the air temperature will adjust 

so that the condensation and evaporation rates are the same. At this condition, the entire 

condensation-evaporation cycle will depend only on the temperature difference between the hot 

and cold plates, irrespective of the temperature distribution of air. 

Figures 3-15 and 3-16 show the temperature comparison at thermocouple locations T3, T4, and 

T5 as indicated in Figure 3-3 for sets 1, 3, and 4. Only three sets of data are presented for 

clarity. The majority of the temperature data have been underpredicted by the nonequilibrium 

model, whereas the match between the experimental and computed data is better for the 

equilibrium model. As with the case of condensation rate, the better match using the equilibrium 

model can be attributed to including volumetric condensation, which is not considered in the 

nonequilibrium model. Consequently, any thermal effect due to exchange of latent heat is not 

accounted for in the nonequilibrium model. In reality, condensate water releases the latent heat 

in the continuous flow resulting in a higher temperature forthe gas phase. This phenomenon is 

accounted for in the equilibrium model, where the latent heat of condensation is absorbed by 

the gas mixture resulting in a higher temperature that is closer to the experimental observation. 

3.3 Summary 

This chapter discussed a combined experimental and computational study to understand the 

condensation-evaporation cycle in an enclosed space. The primary objective of the study was 

to verify and benchmark the numerical technique described in Chapter 2 for further use in 

simulating multiphase natural convection dominated cold-trap process. The experimental study 

explore the effect of volumetric condensation. The results presented showed the velocity, 

temperature, and relative humidity distribution of the system. It was observed that the 
3-12 

Figure 3-15. Comparison of Experimental and Computed 

Values of Temperature Using Nonequilibrium Model 

(Without Volumetric Condensation) [°F = 1.8 × °C + 32]

Figure 3-16. Comparison of Experimental and Computed Values of 

Temperature Using Equilibrium Model 

(With Volumetric Condensation) [°F = 1.8 × T °C + 32]

was conducted in a closely controlled environment, where a condensation-evaporation cycle 

was created within an enclosure using a water source and cold plate at constant temperatures. 

The experimental study was conducted for a range of temperature differences between the cold 

plate and hot water source. Subsequently, a numerical study was conducted to simulate the 

experimental setup. The numerical methodology described in Chapter 2 was used to replicate 

the experiment. Both nonequilibrium and equilibrium multiphase models were studied to 
3-13 

equilibrium model predicted moisture content level bounded by saturation limit unlike the 

nonequilibrium model that predicted supersaturation. Numerical results of the condensate rate 

using both models matched the general experimental trends and patterns, but the equilibrium 

model provided a better prediction due to the availability of additional water vapor at the cold 

wall. Similarly, the numerical simulations provided temperature distributions that compared well 

with experimental observations. However, the equilibrium model provided a better match as 

compared to the nonequilibium model. The temperature values computed using the 

nonequilibrium model were lower as it ignored the latent heat absorption due to volumetric 

condensation process. 
4-1 

4 20-PERCENT SCALE MODEL OF A DRIFT 

This chapter describes an experimental and numerical study conducted using a 20-percent 

drift-scale configuration. The objective of the experimental study was to gain insights into the 

cold-trap process and provide measured data for benchmarking and validating the 

computational fluid dynamics (CFD) simulations. The data collected from the experiments in the 

20-percent drift-scale also provided boundary and initial conditions necessary to perform the 

numerical analysis. Details of the experimental setup and experimental results are documented 

in Manepally, et al. (2006) and Green (2010). 

There are a number of significant differences between the experimental and 

benchmarking exercise described in Chapter 3 and the study described in this chapter. 

The condensation-evaporation cell experiment described in the previous chapter was conducted 

in a closely controlled environment, with minimum heat loss to the environment and stringent 

monitoring of the water source and cold plate temperatures. The condensation cell was also 

configured such that a single two-dimensional convection cell in the vertical plane would form 

within the enclosure. There was little cross flow or heat exchange in the horizontal plane. 

Hence a two-dimensional numerical simulation along the central vertical plane was sufficient to 

model the problem. On the other hand, the 20-percent drift-scale experiment was conducted in 

a relatively open environment to replicate the actual atmosphere of an underground drift. The 

experimental set-up contained a closed thick cylinder, but no special effort was made to restrict 

heat and mass transfer between the experimental system and the surroundings. The flow and 

heat transfer were also expected to be strongly three-dimensional with cross flow convection 

and heat transfer. Another significant difference between the convection cell experiment and 

the 20-percent drift-scale study lies in the mode of heat transfer. Due to relatively low 

temperature difference between the hot water and cold plate, radiation heat transfer was 

not significant in the previous study as described in Section 3. However, in the 20-percent 

drift-scale experiment, radiation heat transfer plays an important role. Overall, the focus of the 

benchmarking study in this chapter is to assess the capability of the CFD model in simulating 

the flow and thermal environment, including the cold-trap process in an emplacement drift 

4.1 The 20-Percent Drift-Scale Experiment 

The experimental setup design and results are discussed in this section. Two distinct sets of 

experiments were conducted. In the first set, the temperature variation due to natural 

convection was studied without any moisture transport. The second set included the effect of 

moisture transport and measured the temperature and relative humidity in the domain. 

4.1.1 Design of the Experiment 

The experimental setup consists of a polyvinylchloride pipe (PVC) closed on each end to 

simulate the enclosed environment and four aluminum analog waste packages as shown in 

Figure 4-1. The pipe that simulates the drift walls is approximately 6.51 m [21.35 ft] in length, 

with internal and external diameters of approximately 1.06 m [3.48 ft] and 1.134 m [3.72 ft], 

respectively. The end caps on the pipe are made of 1.3 cm [0.5 in] thick low thermal 

conductivity Lexan® to minimize heat loss and allow visual observation of the experiment. 
4-2 

Figure 4-1. Photographs of the 20-Percent Scale Drift Natural Convection and 

Cold-Trap Laboratory Model With (a) Polyvinylchloride Pipe, (b) Four Analog 

Waste Packages and Stands, and (c) Waste Packages Inside the Pipe 

With Thermocouples to Measure Surface Temperature 

Figure 4-2 shows the 20-percent drift-scale pipe geometry and dimensions. The axial cross 

section of the pipe shows the offset of the waste packages placed eccentrically into the drift pipe 

and separated by a constant distance of 5 cm [2 in] from the left end of the pipe (hot end wall) 

and between adjacent waste packages. The extreme ends of the pipe were sealed to minimize 

the flow of air in and out of the pipe. The end cap of the pipe closest to the waste packages will 

be referred to as the hot end wall throughout this report; likewise, the other end cap which is far 

away from the waste packages is called the cold end wall. 

(a) 

(b) 

(c) 
4-3 

Figure 4-2. Schematic Drawing of the 20-Percent Drift-Scale Laboratory Model. 

Longitudinal Cross Section Shows the Length of the Tube and Waste Packages and 

Their Location Relative to the Hot End Wall. Four 20-Percent Scale 

Waste Packages (WPA, WPB, WPC, and WPD) Are Placed Eccentrically in the 

Z-Direction. Dimensions Are in Meters. [1 m = 3.28 ft]

Four waste packages are used to simulate the effects of uniform and non uniform heat loads on 

the evolution of the cold-trap process. Each waste package is approximately 30.5 cm [12 in] in 

diameter and approximately 1 m [3.28 ft] in length. Extending internally from one face of each 

waste package is a heating rod. During the experiments, the rod is heated in the presence of a 

vacuum pressure of 689.5 Pa [0.1 psia] inside the analog waste package so that radiation is the 

dominant heat transfer process between the rod and the walls resulting in uniform heating of the 

waste packages. 

For tests with moisture, water was introduced inside the PVC through tubing located right above 

the waste package closest to the hot end (Figure 4-2). Paper towels were located on the top of 

the first waste package (WPA) to (i) distribute the liquid on the surface of the waste package, 

(ii) avoid water dripping off the waste package, and (iii) enable observation as the experiment 

progressed (Figure 4-3). 

4.1.2 System Monitoring 

During the experiments, the temperature was monitored using 120 nested calibrated 

thermocouples on the surface of the waste packages, suspended in the air, and on the inner 

and outer walls of the pipe. These thermocouples were placed in different axial locations that 

are shown in Figure 4-4(a). Each cross section location was instrumented to obtain a 

comprehensive map of the temperature field. An example of the location of thermocouples in a 

cross section of the pipe is shown in Figure 4-4(b). Detailed information on the location of the 

thermocouples can also be found in Manepally, et al. (2006). 

In the following sections of this report, the term “cross sectional average temperature” will refer 

to the average of the individual measured or simulated air temperature values at the different 

locations in a cross section as shown in Figure 4-4a. The number of points in each cross 

section was based on preliminary CFD simulations. Sensors were located in areas where 

significant fluid (air) movement and/or temperature gradients were anticipated. 

6.509 4.200 3.200

3.150

Ø1.060

Ø1.134

2.150

2.100 1.100

1.050

0.050

Package

Waste

Drift

Pipe

WPA WPB WPC WPD

0.013 0.013

Z

X

0.183

0.195

Ø0.305

Cold End Wall Hot End Wall
4-4 

Figure 4-3. Photograph of Moist Area on Top of Waste Package A 

During Experiments With Moisture 

During the testing with moisture, temperature was monitored using thermocouples and 

thermistors. Relative humidity was monitored using a series of sensors at 15 different locations 

along the centerline of the drift and at different heights. The humidity measurements were 

made with commercially available humidity sensors. 

4.1.3 Experimental Conditions 

The thermal and air flow field within the experimental setup were expected to change with time. 

The use of constant boundary conditions and heat rates, however, allowed the system to reach 

a point where average in-drift parameter change was negligible as compared to the initial 

transient response at the beginning of the test. Hence, all of the experimental measurements 

were taken after the initial transient period was over and the system had reached a steady state. 

As mentioned previously, two main scenarios were considered in the testing. For the first 

scenario, the drift was instrumented as described in Section 4.1.2, and cap ends were installed 

at each end of the pipe isolating the ambient air within the pipe from the rest of the room. 

However, as no water was introduced internally, ‘dry’ is used here to identify any experimental 

or numerical results obtained under this first scenario. Only temperatures were measured for 

the ‘dry’ cases. 

The second scenario involved adding water inside the 20-percent drift-scale pipe to represent a 

potential source of liquid in the drift. Liquid water was added on top of the waste package 

closest to the hot end wall (Figure 4-3). During testing, temperature and relative humidity 
4-5 

(a) 

(b) 

Figure 4-4. (a) Schematic Showing the Location of Axial Planes Selected for 

Placing Thermocouples and (b) Example Location of Thermocouples at a Cross 

Section and Calculation of the Cross Sectional Average Air Temperature

were measured at different locations. This part of the analysis aids in understanding both 

temperature and moisture redistribution along the 20-percent drift-scale model. The results 

and analyses related for the second scenario are identified by ‘moist’ in the following sections 

of this report. 

Four different combinations of heating rates applied to the waste packages were selected to 

investigate the effect of the total heat load and distribution on the cold-trap process inside the 

drift pipe (Table 4-1). The nomenclature in Table 4-1 is used in the following sections of this 

report to identify the different configurations of heat load distribution applied during the testing 

and CFD analysis. The four different heat load distributions shown in Table 4-1 were used 

during the tests with dry air, and only two cases, one with uniform {50 W [170.6 BTU/h]} and one 

with nonuniform heat distributions {75-25-25-75 W [266-85-85-266 BTU/h]}, were employed 

during the tests with moisture. 

2.650

5.830

5.130

4.430

3.950

3.700

2.125

1.850

1.600

0.550

A1 BM1 A2 A3 A4 BM2 A5 B1 CM1 C1

WPA WPC WPD WPB

T 5

WP

T 6

T 7

T 2

T 3

T 8

T 9

T 10

8

T2 T3 T5 T6 T7 T8 T9 T10 Cross sectional average+ + + + + + + − =

T 5

WP

T 6

T 7

T 2

T 3

T 8

T 9

T 10

8

T2 T3 T5 T6 T7 T8 T9 T10 Cross sectional average+ + + + + + + − =
4-6 

Table 4-1. Thermal Loads for the 20-Percent Scale Experiments 

Waste 

Package WPA WPB WPC WPD

Total Heat 

Rate (W)* Nomenclature 

In-Drift Fluid Condition: 

Moist (M)H and Dry (D)I

Uniform 

Heating Rate 

(W)* 

50 50 50 50 200 50 W M and D 

75 75 75 75 300 75 W D 

Non-uniform 

Heating Rate 

(W)* 

75 25 25 75 200 75-25-25-75 W M and D 

80 60 40 20 200 80-60-40-20 W D 

*1 W = 3.4 BTU/h 

HM indicates that both temperature and relative humidity were measured during testing and the User Defined 

Functions for volumetric condensation, wall condensation and evaporation in ANSYS-FLUENT®

 were used in the 

numerical simulations. 

ID also implies that only temperatures were measured during testing and the baseline ANSYS-FLUENT®

 solver was 

used without any User Defined Functions. 

4.1.4 Experimental Results 

The majority of the experimental results are discussed along with numerical data in 

Sections 4.2.2.1 and 4.2.2.2. Some highlights of the experimental data are shown in 

Figure 4-5. A plot of the cross sectional average temperatures as a function of the distance 

along the drift is shown in Figure 4-5 for different variants of thermal loads. Average 

temperatures are always higher for the uniform heating case 75 W [266 BTU/h] as the total heat 

load {300 W [1,024 BTU/h]} is higher than for the three other cases {200 W [682 BTU/h]}. All 

the temperature curves resulting from the 200 W [682 BTU/h] thermal load converge towards a 

single value toward the cold end of the drift. However, the 75 W [266 BTU/h] case shows a 

higher temperature value at drift-end. There is, however, substantial variation in the 

temperature pattern over the waste packages due to the difference in thermal load distribution, 

although the total heat load is the same. In general, there is a decrease in temperature from the 

hot end toward the cold end on the order of 2–4 °C [3.6–7.2 °F]. 

4.2 The 20-Percent Drift-Scale Modeling 

The numerical model developed using ANSYS-FLUENT (ANSYS, Inc., 2009a) for simulating the 

20-percent drift-scale experiment is described in this section. Temperature, relative humidity, 

and velocity contours computed from the model are presented at different cross sectional 

locations. Numerical results are also compared with experimental data for average temperature 

and relative humidity values at specific locations. 

4.2.1 Model Development 

A detailed three dimensional model was developed for the 20-percent drift-scale experimental 

setup using ANSYS-FLUENT. The immediate objective of this exercise was to benchmark and 
4-7 

Figure 4-5. Observed Average Cross Sectional Air-Temperature Profile for 

All the Dry Cases With Different Heat Load Distributions Between Individual 

Waste Packages [1 m = 3.28 ft; °F = 1.8 x °C + 32; 1 W = 3.4 BTU/h] 

verify the developed User Defined Functions (UDFs) for modeling condensation. The broader 

objective was to assess the general numerical approach for its intended use in simulating the 

cold-trap process. 

The computational domain used in the model is shown in Figure 4-6. Some of the dimensions 

and boundary conditions are also highlighted in Figure 4-7. The outer boundary of the domain 

is a cylindrical surface of 0.53 m [1.738 ft] diameter and a thickness of 0.023 m [0.007 ft], that 

represents the drift. A numerical study was conducted to see the effect of the outer wall 

boundary condition on the temperature distribution pattern. Two different boundary conditions 

were tested: (i) a specified temperature at the outermost surface and specified wall thickness, 

where the PVC pipe shell was not included in the computation and (ii) a convective flow 

boundary condition at the outermost surface with a specified ambient temperature, where the 

pipe thickness was included in the computation. It was observed that a convection boundary 

condition on the outermost surface of the PVC pipe provided the best estimate and this 

approach was adopted for rest of the modeling exercises. The ambient temperature used in the 

calculation was obtained from the experimental observations. 

The emplaced 20-percent scale waste packages have a length of 1 m [3.28 ft] and a diameter of 

0.365 m [1.198 ft]. The waste package volumes were included in the computation. It was 

assumed in the computation that the waste packages were solid objects with a uniform thermal 

conductivity and has a uniformly distributed constant volumetric source. The thermal load 

specified in the waste packages was in accordance with the experimental input and was varied 

for different tests to reflect the loading used in the experiment. The air density was calculated 

using a Boussinesq approximation which assumes that density is a function of temperature and 

25

26

27

28

29

30

31

32

33

34

35

01234567

Distance from Hot End Wall of Drift (m)

Average Temperature (°C)

50w

75w

75-25-25-75

80-60-40-20

WPA WPD WPB WPC
4-8 

Figure 4-6. Schematic of the Three-Dimensional ANSYS-FLUENT® Model 

Showing Relevant Dimensions [1 m = 3.28 ft]

Figure 4-7. Schematic of the Three-Dimensional ANSYS-FLUENT® Model Showing 

Dimensions and Boundary Conditions [1 m = 3.28 ft]

is independent of pressure. This assumption is valid as the flow velocity was expected to be 

low and variation in pressure was expected to be minimal. The ‘dry’ and ‘moist’ air properties 

used in the calculation werw obtained from Incropera and DeWitt (1996) and Reynolds (1979). 

The mesh elements used in the model are shown in Figure 4-8. Mostly hexahedral elements 

are used to increase computational efficiency, but some tetrahedral elements were needed to 

effectively mesh small gaps. Figure 4-9 shows the mesh elements used for the waste packages 

and end caps. It also highlights the cross sectional distribution of the mesh. A weep algorithm 

was used to generate the volumetric mesh from a surface distribution. As shown, a finer mesh 

resolution was used near the solid walls to resolve hydrodynamic and thermal boundary layers. 

PVC Pipe Length 6.5344 m 

Waste Package Length 1.0 m

WPD WPC WPB WPA

Diameter 0.365 m 

Water Source 

Top Part of 

WPA 

PVC Pipe 

Diameter 

0.53 m 

PVC Pipe Thickness 0.023m

End Cap 

Thickness 0.023 m 

Specified Volumetric Heat 

Source in Waste Packages 

WPA WPB WPC and WPD 

Convection Heat 

Transfer on Outer 

Boundaries 

WPA 

WPB

WPC

WPD
4-9 

Figure 4-8. Volumetric Mesh for the Computing Domain Used 

in the ANSYS-FLUENT® Model

Figure 4-9. Highlight of Meshing Over Waste Packages and End Caps 

Used in the ANSYS-FLUENT® Model
4-10 

4.2.2 Model Results 

Modeling results are presented in two sections. Section 4.2.2.1 presents results for the ‘dry’ 

cases, where no liquid water was available for evaporation. All four scenarios for the thermal 

loading described in Table 4-1, which includes two uniform {50 and 75 W [170.6 and 

266 BTU/h]} and two nonuniform {75-25-25-75 and 80-60-40-20 W [266-85-85-266 and 

273-205-137-68 BTU/h]} heat load distributions, were simulated. For the ‘dry’ cases, all three 

modes of heat transfer, conduction, convection and radiation were included in the model. As 

mentioned in Section 2.3, the discrete ordinate model was used to calculate radiation. 

Temperature contours and velocity vectors at different cross sectional locations are presented 

for the ‘dry’ case with a 50 W [170.6 BTU/hr] thermal load. In addition, cross section averaged 

temperature values from the computed data are compared with available experimental data for 

all the simulated runs. 

Section 4.2.2.2 presents results for the ‘moist’ case where a source of liquid water was 

specified to exist on the upper half of the waste package WPA outer surface and was 

available for evaporation. As indicated in Table 4-1, two different thermal loading scenarios 

were considered for the ‘moist’ case; one uniform loading of 50 W [170.6 BTU/hr] and one 

nonuniform loading of 75-25-25-75 W [266-85-85-266 BTU/hr]. For the ‘moist’ cases, relative 

humidity and temperature distributions at different cross sectional locations are presented with a 

50 W [170.6 BTU/hr] thermal load. Also, the section averaged temperature and relative 

humidity data obtained from the simulation are compared with the experimental results for all 

the simulated runs. 

4.2.2.1 Simulation of Dry Cases 

Figures 4-10(a) and (b) show the streamlines obtained from the 50 W [170.6 BTU/hr] ‘dry’ 

simulation case. The streamline ribbons track the flow path and the colors indicate velocity 

magnitude. The streamlines highlight the general velocity pattern generated due to buoyancy 

driven natural convection. In general, a large axial circulation loop can be observed from both 

Figures 4-10(a) and (b). In this loop, air moves from the hot to the cold end of the drift in the 

space above the waste packages. Air convects back from the cold to the hot end through a gap 

between waste packages and the bottom wall. Streamlines in Figure 4-10(a), however, also 

indicate the three dimensionality of the flow field, due to the presence of strong vertical currents 

in planes between waste packages. Away from the waste packages, towards the cold end of 

the tube, the flow is mostly in the axial direction and exhibits a two dimensional pattern that is 

mostly highlighted in Figure 4-10(b). 

Axial velocity magnitude contours and velocity vectors at different cross sectional planes will 

help to further explore the pattern of the velocity distribution in different regions. Figure 4-11 is 

a schematic that shows different cross sectional planes, where velocity and temperature 

contours are presented. 

Figures 4-12 through 4-16 show the velocity distribution in the cross sections defined in Figure 

4-11 for four different thermal loading conditions. The contours in these diagrams represent 

axial velocity magnitude. Figure 4-12 shows the results for the 50 W [170.6 BTU/hr] uniform 

thermal loading scenario and will be discussed here as the baseline result. It is clear from the 

vector plots of axial planes 1 through 4 that due to the thermal effect, strong velocity plumes are 

created around the waste packages that convect in the vertical direction. This finding agrees 

with the experimental observations of Kuehn and Goldstein (1976 and 1978), where they 
4-11 

Figure 4-10. Airflow Streamlines for a Uniform 50 W [170.6 BTU/hr] 

Heat Load [1m/s=3.28 ft/s]

observed a similar flow pattern. The velocity distribution pattern in the vertical YZ plane 

confirms the existence of these vertical convection currents and also highlights the two 

dimensional flow pattern away from the waste packages towards the cold end of the tube. A 

high axial velocity region near the top surface can be observed in the vertical YZ plane after the 

flow leaves the waste package WPD. The combined air mass that convected vertically upwards 

due to the thermal effect of the waste packages, tends to move axially once it leaves the waste 

package area and creates a high velocity region. 

Figures 4-13 through 4-15 highlight the effect of uniform and nonuniform thermal loading on the 

velocity distribution pattern. As is clear from Figure 4-13, the high load waste packages {75 W 

[266 BTU/hr] for WPA and WPD} creates a stronger vertical movement evident from the higher 

velocity magnitude compared to the low load waste packages {25 W [85 BTU/hr]} for WPB and 

WPD). Similar observations can be made in Figure 4-15, where the largest vertical velocities 

are observed above WPA that has the maximum heat load of 80 W [275 BTU/hr] and the 

(a) 

(b) 

Hot End 

Cold End 

WPD WPC WPB WPA 
4-12 

Figure 4-11. Schematic of Cross Sectional Planes 

Figure 4-12. Velocity Vectors and Axial Velocity Contours for the Dry Simulation 

With a Thermal Load of 50 W [170.6 BTU/hr; 1m/s=3.28 ft/s]

Hot End 

Cold End 

Vertical YZ Plane 

Axial Plane-1 along Mid 

Axial Length of WPA 

Axial Plane-2 Along the 

Middle of the Gap between 

WPA and WPB 

Axial Plane-3 along Mid Axial 

Length of WPD 

Axial Plane-4 Middle of Cold and End of WPD 

WPC

WPB

WPD

WPA

Axial Plane-1 

Axial Plane-2 

Axial Plane-3 

Axial Plane-4 

Vertical YZ Plane 

WPD 

WPC

WPB WPA 
4-13 

Figure 4-13. Velocity Vectors and Axial Velocity Contours for the Dry Simulation 

With Thermal Loads of 75-25-25-75 W [266-85-85-266 BTU/hr; 1m/s=3.28 ft/s] 

Figure 4-14. Velocity Vectors and Axial Velocity Contours for the Dry Simulation 

With a Thermal Load of 75 W [266 BTU/hr; 1m/s=3.28 ft/s]

Axial Plane-1 

Axial Plane-2 

Axial Plane-3 

Axial Plane-4 

Vertical YZ Plane 

WPD 

WPC

WPB WPA 

Axial Plane-1 

Axial Plane-2 

Axial Plane-3 

Axial Plane-4 

Vertical YZ Plane 

WPD 

WPC

WPB

WPA 
4-14 

Figure 4-15. Velocity Vectors and Axial Velocity Contours for the Dry Simulation 

With Thermal Loads of 80-60-40-20 W [273-205-137-68 BTU/h; 1m/s=3.28 ft/s]

velocity magnitude gradually decreases to a minimum above WPD that has the minimum heat 

load of 20 W [68 BTU/hr]. On the other hand Figure 4-14 shows that the vertical convection 

generated due to waste package heating has similar magnitude as all the waste packages have 

the same thermal loading of 75 W [266 BTU/hr]. Figures 4-13 through 4-15 also show that the 

relative location and magnitude of the high axial velocity region mentioned previously also 

depends on the thermal load distribution. The largest axial high velocity zone with maximum 

velocity magnitude was generated by the most uneven thermal load distribution; 80-60-40-20 W 

[273-205-137-68 BTU/h] as shown in Figure 4-15. The uniform distribution case with minimum 

heat load 50 W [170.6 BTU/hr], Figure 4-12 generated the smallest high velocity zone. 

Figure 4-16 shows the temperature calculated at different cross sectional locations illustrated in 

Figure 4-12. The temperature distribution in the axial planes for all four ‘dry’ cases exhibited the 

same pattern of temperature distribution, though the local values were dependent on the 

thermal loading. A comparison of temperature contours in the YZ plane among Figures 4-16(a) 

through (d) shows that the temperature values were dependent on the thermal loading of the 

adjacent waste package. For example Figure 4-16(d) shows that the temperature above WPA 

which has a load of 80 W [273 BTU/hr] is higher compared to temperature above WPD, that has 

a load of 20 W [68BTU/hr]. This phenomenon illustrates the fact that the primary heat transfer 

mode in the vicinity of the waste package is two-dimensional vertical buoyancy driven 

convection and as illustrated previously in Figures 4-12 through 4-15, the convective velocity 

magnitude is dependent on the thermal loading. A qualitative comparison of temperature 

contours in the YZ plane between different cases also illustrates that axial velocity has only a 

Axial Plane-1 

Axial Plane-2 

Axial Plane-3 

Axial Plane-4 

Vertical YZ Plane 

WPD 

WPC

WPB

WPA 
4-15 

Figure 4-16. Temperature Contours Dry Simulations [°F = 1.8 × K–459.4] 

minor impact on heat transfer from waste packages and any axial convection is likely to be 

small compared to vertical convection. 

Waste package temperatures calculated for the ‘dry’ cases are illustrated in Figure 4-17. As is 

evident from Figure 4-17, high thermal load waste packages exhibit higher surface 

temperatures compared to the low thermal load waste packages. In addition, a high thermal 

load waste package had minimal thermal impact on an adjacent low thermal load waste 

package. Hence there has been little or no heat exchange between waste packages due to 

radiation and it can be concluded that waste package temperatures are mostly a function of the 

thermal load it contains. 

Axial Plane-1

Vertical YZ Plane 

Axial Plane-2 

Axial Plane-3

Axial Plane-4 

(a) 50 W [170.6 BTU/hr] (b) 75-25-25-75 W [266-85-85-266 BTU/hr] 

(c) 75 W [266 BTU/hr] (d) 80-60-40-20 W [273-205-137-68 BTU/hr]

WPB WPA

WPC 

WPD 
4-16 

Figure 4-17. Temperature Distribution on Waste Package Surfaces 

[1 m =3.28 ft; °F = 1.8 × K–459.4]

Figures 4-18(a) through (d) show comparisons of experimental temperature data with computed 

results. These figures illustrate the section averaged temperature data as a function of the 

distance along the drift for all four ‘dry’ cases. Location of the axial cross sections are shown in 

Figure 4-4(a) and the section average temperature calculation procedure is explained in Figure 

4-4(b). Details on the instrumentation in each section are reported in Manepally, et al. (2006) 

and Green (2010). The temperature results predicted here include all modes of heat transfer, 

including radiation, but as the domain does not have any water source, moisture transport was 

excluded. In general, it can be seen from Figure 4-18 that the computations were able to 

replicate the general trend and pattern of the experimentally observed temperature variation. 

Most of the temperature peaks observed in the experiment were captured in the computations. 

This indicates that the computations were successful in obtaining the spatial temperature 

variation within the domain. But it can also be observed that for all the cases, the computed 

results consistently overpredict the average temperature values. The difference between the 

computed and experimental data, however, is not constant and depends on the point of 

(a) 50 W [170.6 BTU/hr] (b) 75-25-25-75 W [266-85-85-266 BTU/hr] 

(c) 75 W [266 BTU/hr] (d) 80-60-40-20 W [273-205-137-68 BTU/hr] 

WPA 

WPB 

WPC 

WPD 
4-17 

Figure 4-18. Cross Section Averaged Temperature Distribution Along the 

Axial Distance [1 m = 3.28 ft; °F = 1.8 × °C + 32] 

observation and ranges between 2 °C [3.6 °F] to 0.3 °C [0.54 °F]. The possible reason for the 

discrepancy between the calculated and experimental data is discussed below in connection 

with Figure 4-19. 

Figure 4-19 shows a comparison between computed and experimentally obtained averaged 

waste package temperature distributions. Contrary to what has been observed in Figure 4-18, 

the computed data underpredicts the average waste package temperature. The difference 

between computed and experimental values ranges between 2.3 °C [4.14 °F] to 0.4 °C 

[0.72 °F]. Results from Figures 4-18 and 4-19 show that in the computation a model excess 

heat is released from the waste packages and warms up the adjacent air, which causes 

underprediction of waste package temperatures and overprediction of section averaged air 

(a) 50 W [170.6 BTU/hr] (b) 75-25-25-75 W [266-85-85-266 BTU/hr]

(c) 75 W [266 BTU/hr] (d) 80-60-40-20 W [273-205-137-68 BTU/hr] 
4-18 

 

 

Figure 4-19. Temperature Distribution on Waste Package Surfaces 

[1 m = 3.28 ft; °F = 1.8 × °C + 32]

 

temperatures. This result can also be attributed to uncertainty in material property specification 

such as emissivity that controls the heat release rate from the waste packages. 

To have an overall assessment of the model, temperatures obtained from the numerical solution 

are plotted against experimental data in Figure 4-20. It is clear that for the majority of the cases, 

the computed results overpredict temperature. The deviation in results, however, does not 

exceed 3 °C [5.4 °F]. 

4.2.2.2 Simulation of Moist Cases 

This section discusses the computational results for the ‘moist’ case simulations as identified 

in Table 4-1. Experimental investigation for the ‘moist’ case was performed for one 

uniform thermal loading of 50 W [170.6 BTU/h] and one nonuniform loading of 75-25-25-75 W 

[266-85-85-266 BTU/h]. To get an understanding of the thermal effect on moisture distribution 

(b) 75-25-25-75 W [266-85-85-266 BTU/hr]

(c) 75 W [266 BTU/hr] (d) 80-60-40-20 W [273-205-137-68 BTU/hr] 

(a) 50W [170.6 BTU/hr] 
4-19 

Figure 4-20. Comparison of Experimental and Computed Temperature Data for all Dry 

Cases With Uniform and Nonuniform Heat Load Distributions 

[°F = 1.8 × °C + 32]

inside the drift, the numerical study was done for all four heat loads mentioned in Table 4-1, 

though experimental data was available for only two cases. 

The standard ANSYS-FLUENT solver along with the UDFs developed to model wall 

condensation-evaporation and volumetric condensation was used to obtain results for the 

‘moist’ cases. Both equilibrium and nonequilibium models in the simulation provided very 

similar results for flow and temperature fields. This is because the volumetric condensation 

within the system was almost absent as water vapor concentration mostly remained at or below 

the saturation limit. The amount of water vapor generated due to evaporation of dripping water 

was relatively modest compared to the condensation cell chamber described in Chapter 3. 

Hence, results obtained from nonequilibrium model are presented and discussed in this section. 

Thermocouples and humidity sensors were placed at twelve spatial locations for the ‘moist’ 

cases for both uniform and nonuniform loading. The positions of these sensors are shown 

in Figures 4-21 and 4-22. As shown in these figures, all probes were placed along the 

mid-vertical YZ plane {x=0 m [0 ft]}, and at five axial planes. Locations of these axial planes 

are shown in Figure 4-21. Figure 4-22 highlights the relative position of these sensors in the 

mid-vertical plane. 

Figure 4-23 shows the axial variation of experimental and computed cross section averaged 

temperature for both the uniform 50 W [170.6 BTU/hr] and non uniform 75-25-25-75 W 

[266-85-85-266 BTU/hr] thermal loading ‘moist’ case. The simple sectional averaging method 

explained in Figure 4-4(b) is used for average temperature calculations. It can be observed in 

Figure 4-23 that the computational results predicted the trend of temperature variation, but have 

generally overpredicted the temperature. This was also observed for the ‘dry’ cases and 
4-20 

Figure 4-21. Axial Location of the Relative Humidity Probes for the 

Moist Case Experiments [1 m = 3.28 ft]

Figure 4-22. Relative Position of 12 Thermocouples and Relative 

Humidity Sensors in the Middle Vertical Plane 

highlighted in Figure 4-18. As discussed in Section 4.2.2.1, this can be attributed to 

uncertainties in specifying exact material properties and other modeling parameters. It should 

also be noted that the experimental data was collected for a very limited number of points and 

does not provide a comprehensive map of the temperature variation in the axial direction. 

Deviation of computed data from experimental observations can also be caused by 

measurement uncertainties and localized unsteady flow effects. 

Z= 0.602 m

Z= 1.512 m 

Z= 3.193 m

Z= 3.613m

Z= 5.527m

Point-1 

Point-4 Point-2 

Point-3 Point-5

Point-8 Point-6

Point-9 Point-7

Point-12 

Point-10 

Point-11 
4-21 

Figure 4-23. Cross Section Averaged Temperature Distribution Along the 

Axial Distance for Moist Cases [1 m = 3.28 ft; °F = 1.8 × °C + 32]

Figure 4-24 shows the overall performance of the solver in predicting air temperature at different 

locations within the 20-percent drift-scale with uniform and nonuniform thermal loading for 

‘moist’ cases. It can be seen that the computations generally overpredicted temperature values, 

but for most of the observation points, the deviation did not exceed 3 °C [5.4 °F]. 

Figure 4-24. Comparison of Experimental and Computed Temperature Data for All 

Moist Cases With Uniform and Nonuniform Heat Load Distributions [°F = 1.8 × °C + 32]

(a) 50 W [170.6 BTU/hr] (b) 75-25-25-75 W [266-85-85-266 BTU/hr]
4-22 

Figure 4-25 shows the relative humidity contours for all four thermal loads. The liquid water 

source is on the upper half of waste package WPA. A high humidity area in the surroundings of 

WPA can be observed from the contours in axial-plane-1, that can be attributed to the 

evaporation of liquid water from the surface of WPA. The level of humidity and water vapor 

mass concentration, however, is dependent on the thermal load of WPA. A qualitative 

comparison of 4-25(a) and 4-25(d) shows that a thermal load of 50 W [170.6 BTU/hr] sustains 

more humidity compared to a heat load of 80 W [173 BTU/hr]. This was expected, as the air 

Figure 4-25. Relative Humidity Contours Obtained From Moist Simulations 

Axial Plane-1

Vertical YZ Plane 

Axial Plane-2 

Axial Plane-3

Axial Plane-4 

(a) 50 W [170.6 BTU/hr] (b) 75-25-25-75 W [266-85-85-266 BTU/hr]

(c) 75 W [266 BTU/hr] (d) 80-60-40-20 W [273-205-137-68 BTU/hr] 

WPA WPB 

WPC 

WPD 
4-23 

temperature in the vicinity of a waste package is mostly controlled by the vertical convection 

current generated due to the waste package heat and is not significantly influenced by the axial 

convection or inter-waste package radiation as discussed in Section 4.2.2.1. Hence air around 

a high thermal load waste package has a higher temperature, as illustrated in Figure 4-16, and 

has a higher saturation limit that is manifested as lower relative humidity. Similar dependence 

of relative humidity on the waste package heat load can be observed for the rest of the waste 

packages as well. The relative humidity contours on axial planes show that the air below the 

waste package has higher relative humidity. This can also be attributed to the convection 

pattern of hot air as illustrated in Figures 4-12 through 4-15 that results in a low temperature 

zone below the waste package as highlighted in Figure 4-16. The relative humidity near the 

cold end of the domain is also higher due to relatively cooler temperature in that region. 

Figure 4-26 compares the experimental and computed relative humidity data. The axial 

variation of relative humidity for the ‘moist’ uniform 50 W [170.6 BTU/hr] thermal loading is 

presented in Figure 4-26(a) and for ‘moist’ nonuniform 75-25-25-75 W [266-85-85-266 BTU/hr] 

thermal loading is presented in Figure 4-26(b). In general, the relative humidity initially 

decreases with axial distance from the hot end and subsequently increases with axial 

distance away from the waste packages. Deviation of this pattern for the 75-25-25-75 W 

[266-85-85-266 BTU/hr] ‘moist’ case at z = 1.512 m [4.96 ft] was caused by a malfunctioning 

relative humidity probe and should be ignored. For both uniform and nonuniform thermal 

loading, computed results overpredict relative humidity near the hot end and then underpredict 

data away from the hot end. The overprediction of relative humidity near the hot end can be 

associated with uncertainty in estimation of liquid water mass. The evaporation model uses 

surface area and surface temperature as input to calculate the steady vaporization rate. It does 

not take into account the effect of limited liquid water dripping that was done in the experiment 

 

Figure 4-26. Cross Section Averaged Relative Humidity Distribution 

Along the Axial Distance for Moist Cases [1 m = 3.28 ft] 

(a) 50 W [170.6 BTU/hr] (b) 75-25-25-75 W [266-85-85-266 BTU/hr] 
4-24 

or any unsteady effect associated with it. As a result, the computational model allows more 

water to vaporize and causes higher relative humidity near the hot end. Similarly, the difference 

in condensation rate between the experimental observation and computational estimate is 

responsible for the deviation in relative humidity away from the hot end. 

Figure 4-27 shows the overall variation of the computed relative humidity data with the 

experimental data and it is shown that the calculated relative humidity values lie within 

10 percent of the experimental data except at a single location. As mentioned previously, 

this deviation was caused by instrumentation malfunction and should be neglected. 

4.3 Summary 

A 20-percent scale model of the proposed emplacement drift at Yucca Mountain was developed 

to study the heat and moisture transportation due to radioactive decay heat and liquid water. 

The experiment was also conducted to benchmark and test the customized UDFs developed for 

modeling the condensation and evaporation processes to be used in conjunction with the 

baseline ANSYS-FLUENT CFD solver. 

Two sets of tests were performed to understand the waste package heat driven natural 

convection and moisture transport separately. In the first set of tests, no moisture was allowed 

inside the domain and tests were conducted for two uniform and two nonuniform thermal 

loading scenarios. In the second set of tests, liquid water was dripped on the first waste 

package which evaporated to form moist air within the system. One uniform and one 

nonuniform thermal loading was considered for the ‘moist’ system. 

Figure 4-27. Comparison of Experimental and Computed Relative Humidity for 

the Moist Cases With Uniform 50 W [170.6 BTU/hr] and Nonuniform 75-25-25-75 W 

[266-85-85-266 BTU/hr] Thermal Load Distributions 
4-25 

A three dimensional model was developed using ANSYS-FLUENT for the 20-percent scaled 

model. For the ‘dry’ cases, the baseline solver was used to model conduction, natural 

convection and thermal radiation heat transfer. The discrete ordinate radiation model was used 

for radiation modeling. For the ‘moist’ cases, wall condensation and evaporation was simulated 

using the customized UDFs. As some initial simulations exhibited negligible volumetric 

condensation, only the nonequilibium model was used for the ‘moist’ case simulations. 

The flow field obtained from the ‘dry’ cases showed a strong two-dimensional convection cell 

in the axial plane around a waste package, where a hot air plume moved vertically upwards. 

The velocity magnitude of air in these localized cells was directly dependent on the waste 

package thermal load and was not significantly dependent on the heat load of neighboring 

waste packages. Away from the heated waste packages, the two dimensional axial convection 

current was dominant. The flow below the waste packages is not representative of that of the 

20-percent drift-scale model because the metallic stands of the waste packages were not 

included in the CFD model. In the absence of the stands, the flow below the waste packages is 

mainly axial, moving from the cold end to the hot end. A comparison of experimental and 

computed data showed that ANSYS-FLUENT generally overpredicted the temperatures, but the 

deviation did not exceed 3 °C [5.4 °F] 

The ‘moist’ case simulations showed that the relative humidity and the water vapor 

concentration around the waste package were mostly dependent on the heat load of the 

waste package and showed minimal effect of axial convection current or the heat load of 

neighboring waste packages. The computed solution overpredicted relative humidity near 

the water source at the hot end and underpredicted relative humidity away from the waste 

packages near the cold end. The temperature values obtained from the simulations were 

consistently higher compared to the experimental data, but the maximum difference was in the 

range of 3 °C [5.4 °F]. 
5-1 

5 SUMMARY 

The in-drift repository environment is a result of a complex interactions among a number of 

thermal, hydrological and chemical processes. The heat transfer mechanism inside the drift 

involves conduction, natural convection, and radiation as well as latent heat transfer. Liquid 

water inside the drift evaporates from hotter surfaces, convects to a relatively cooler location 

and condenses to form liquid water. This redistribution of water can affect the long-term 

integrity of the engineered barrier system. 

General computational fluid dynamics (CFD) tools have been extensively used to simulate flow 

and thermal fields for a variety of applications. The broad objective of the study was to develop 

a general CFD modeling framework that is useful for simulating the flow, heat and mass transfer 

within an emplacement drift. The modeling exercise was supported by experimental studies 

that provided valuable data for validation and benchmarking. 

A general purpose CFD tool ANSYS-FLUENT® has been adapted for simulation of in-drift flow, 

heat and mass transfer processes. The standard solver was used to simulate the natural 

convection and multimode heat transfer processes including radiation. The solver was also 

capable of tracking water vapor movement within the domain by solving the species transport 

equations. Customized tools for modeling surface condensation, volumetric condensation and 

surface evaporation were developed and linked with the solver as User Defined Functions 

(UDFs). Two different approaches to model volumetric condensation using equilibrium and 

nonequilibium models were developed and tested. Data from supporting experimental studies 

were used for verification and benchmarking of computed data. 

Two different sets of combined experimental and computational investigations were described 

in this report. The first study was conducted in a closely controlled environment, where a 

condensation-evaporation cycle was created within an enclosure using a water source and 

cold plate at constant temperatures. The experimental study was conducted for a range of 

temperature differences between the cold plate and hot water source. Subsequently, a 

numerical study was conducted to simulate the experimental setup. The results showed the 

velocity, temperature, and relative humidity distribution of the system. It was observed that 

the equilibrium model provided a better match with the experimental condensation rate 

and temperatures. 

The second experimental and computational study was conducted using a 20-percent scale 

model of the proposed emplacement drift at Yucca Mountain. Two different sets of experiments 

were conducted. In the first set, no liquid water was introduced into the system and the 

experiment was conducted for a uniform and nonuniform thermal loading of waste packages. 

Temperature measurements were taken for a number of spatial locations within the 

experimental setup. A numerical model was developed to simulate the experimental setup with 

no moisture using ANSYS-FLUENT. Computed velocity patterns showed existence of strong 

two-dimensional convection cells in the vicinity of the waste packages, where the air moves 

vertically upwards. This motion was found to be controlled by the waste package heat load and 

was almost independent of the axial convection or the thermal loading of any neighboring waste 

packages. A comparison of experimental and computed data showed that ANSYS-FLUENT 

generally overpredicted the temperatures, but the deviation did not exceed 3 °C [5.4 °F]. The 

second set of experimental and computational studies included the effect of liquid water. A 

source of liquid water on the upper surface was established. The liquid water subsequently 

evaporated and was transported to a cooler location by natural convection. The experimental 
5-2 

study measured the relative humidity and temperature data. Simulation results showed that like 

velocity, relative humidity and water vapor concentration around the waste package was mostly 

dependent on the heat load of the waste package. The computed solution overpredicted 

relative humidity data near the water source at the hot end and underpredicted relative humidity 

away from the waste packages near the cold end. The computed results overpredicted 

temperature values, but the deviation did not exceed 3 °C [5.4 °F]. 

From the comparison of experimental and numerical results from the above mentioned studies, 

it can be concluded that the ANSYS-FLUENT baseline solver was able to predict the general 

trend and pattern of the convective airflow and heat transfer. The deviation in temperature 

predicted between the computed and experimental data was bounded within a small range. 

The relative humidity and condensation rate results obtained from the combined baseline 

solver along with the UDFs provided results comparable to the experiments, but there was 

some deviation of predicted values from experimental data. Overall results indicate the 

modeling framework and supplementary customized UDFs developed using ANSYS-FLUENT 

can be useful in investigating the in-drift repository environment to assess the flow and 

thermal parameters. 
6-1 

6 REFERENCES 

ANSYS, Inc. “FLUENT® User’s Guide.” Version 12.1. Canonsburg, Pennsylvania: 

ANSYS, Inc. 2009a. 

–––––. “ANSYS–WORKBENCH User’s Guide.” Version 12.1. Canonsburg, Pennsylvania: 

ANSYS, Inc. 2009b. 

–––––. “ANSYS–DESIGN MODELERS User’s Guide.” Version 12.1. Canonsburg, 

Pennsylvania: ANSYS, Inc. 2009c. 

–––––. “ANSYS–MESHING User’s Guide.” Version 12.1. Canonsburg, Pennsylvania: 

ANSYS, Inc. 2009d. 

–––––. “ANSYS–CFD–POST User’s Guide.” Version 12.1. Canonsburg, Pennsylvania: 

ANSYS, Inc. 2009e. 

–––––. “FLUENT® Theory Guide.” Version 12.1. Canonsburg, Pennsylvania: ANSYS, Inc. 

2009f. 

–––––. “FLUENT® User Defined Functions Guide.” Version 12.1. Canonsburg, Pennsylvania: 

ANSYS, Inc. 2009g. 

Bechtel SAIC Company, LLC. “In-Drift Natural Convection and Condensation.” 

MDL–EBS–MD–000001. Rev. 00. Las Vegas, Nevada: Bechtel SAIC Company, LLC. 2004. 

Bell, B. “Application Brief: Film Condensation of Water Vapor.” Lebanon, New Hampshire: 

Fluent, Inc. 2003. 

Bird, R.B., W.E. Stewart, and E.N. Lightfoot. Transport Phenomena. New York City, New York: 

John Wiley & Sons. 1960. 

Birkholzer, J.T. and G.L. Danko. “Technical Work Plan for: Thermal-Hydrological Near-Field 

Model Studies and Impact of Natural Convection on Seepage.” OSTI–LBNL–TWP–000006. 

Rev. 00. Las Vegas, Nevada: DOE, Office of Civilian Radioactive Waste Management, Office 

of Science and Technology and International. 2005. 

Brouwers, H.J.H. and A.K. Chester. “Film Models for Transport Phenomena With Fog 

Formation: The Classical Film Model.” International Journal of Heat and Mass Transfer. 

Vol. 35. pp. 1–11. 1992a. 

–––––. “Film Models for Transport Phenomena With Fog Formation: The Fog Film Model.” 

International Journal of Heat and Mass Transfer. Vol. 35. pp. 13–28. 1992b. 

Buscheck, T.A., L.G. Glascoe, and K.H. Lee. “Validation of the Multiscale Thermohydrologic Model 

Used for Analysis of a Proposed Repository at Yucca Mountain.” Journal of Containment Hydrology. 

pp. 62–633. 2003. 

Carey, V.P. Liquid-Vapor Phase-Change Phenomena. 2nd Edition. New York City, New York: 

Hemisphere Publishing Corporation. 2007. 
6-2 

Das, K. and D. Basu. “Validation of Computational Fluid Dynamics Modeling Approach to 

Evaluate VSC-17 Dry Storage Cask Thermal Designs.” San Antonio, Texas: CNWRA. 2007. 

Das, K., D. Basu, and K. Axler. “Validation of Modeling Approach To Evaluate Transportation, 

Aging, and Disposal Canister Thermal Designs.” San Antonio, Texas: CNWRA. 2008. 

Danko, G. and D. Bahrami. “Barometric Pressure Variation Effects on Emplacement Drift 

Environmental Conditions at Yucca Mountain.” Proceedings of the 11th International High-Level 

Radioactive Waste Management Meeting, Las Vegas, Nevada, April 30–May 4, 2006. 

La Grange Park, Illinois: American Nuclear Society. pp. 353–360. 2006. 

–––––. “Heat and Moisture Flow Simulation with Multiflux.” ASME Heat Transfer/Fluids Summer 

Conference, Charlotte, North Carolina. New York City, New York: ASME International. 2004. 

–––––. “Powered, and Natural, Passive Ventilation at Yucca Mountain.” Proceedings of the 

8th International High-Level Radioactive Waste Management Meeting, Las Vegas, Nevada, 

March 30–April 2, 2003. La Grange Park, Illinois: American Nuclear Society. pp. 683–690. 2003. 

Fedors, R.W., S. Green, D. Walter, G. Adams, D. Farrell, and S. Svedeman. “Temperature and 

Relative Humidity Along Heated Drifts With and Without Drift Degradation.” CNWRA 2004-04. 

San Antonio, Texas: CNWRA. 2004. 

Fedors, R.W., D.B. Walter, F.T. Dodge, S.T. Green, J.D. Prikryl, and S.J. Svedeman. 

“Laboratory and Numerical Modeling of the Cold-Trap Process.” San Antonio, Texas: CNWRA. 

2003a. 

Fedors, R.W., G.R. Adams, C. Manepally, and S.T. Green. “Thermal Conductivity, Edge 

Cooling, and Drift Degradation—Abstracted Model Sensitivity Analyses for Yucca Mountain.” 

San Antonio, Texas: CNWRA. 2003b. 

Flow Science, Inc. “FLOW-3D® User’s Manual.” Version 9.0. Santa Fe, New Mexico: 

Flow Science, Inc. 2005. 

Fox, R.J., P.F. Pederson, M.M.L. Corradini, and A.P. Pernsteiner. “Free Convective 

Condensation in a Vertical Enclosure.” Nuclear Engineering and Design. Vol. 177, Nos. 1–3. 

pp. 71–89. 1997. 

Green, S. “Cold Trap CFD Modeling.” Scientific Notebook Number 536E. San Antonio, Texas: 

CNWRA. 2010. 

Green, S. and C. Manepally. “Software Validation Test Plan for FLOW-3D Version 9.” Rev. 1. 

San Antonio, Texas: CNWRA. 2006. 

Green, S.T., D.B. Walter, R.W. Fedors, and F.T. Dodge. “A Model for Moisture Transport in a 

High-level Radioactive Waste Repository Drift.” Proceedings of the Society for Mining, 

Metallurgy, and Exploration Annual Meeting, February 23–25, 2004. Denver, Colorado: Society 

for Mining, Metallurgy, and Exploration. Published on CD ROM. 2004. 

Hijikata K. and Y. Mori. “Forced Convective Heat Transfer of a Gas With Condensing Vapor 

Around a Flat Plate.” Heat Transfer Japanese Research. Vol. 2. pp. 81–101. 1973. 
6-3 

Incropera, F.P. and D.P. Dewitt. Fundamentals of Heat and Mass Transfer. 4th Edition. 

New York City, New York: John Wiley & Sons, Inc. 1996. 

Kang, H.C. and M.H. Kim. “Characteristics of Film Condensation of Supersaturated Steam–Air 

Mixture on a Flat Plate.” International Journal of Multiphase Flow. Vol. 25, No. 8. 

pp. 1,601–1,618. 1999. 

Kuehn, T.H. and R.J. Goldstein. “An Experimental Study of Natural Convection Heat Transfer in 

Concentric and Eccentric Horizontal Cylindrical Annuli.” ASME Journal of Heat Transfer. 

Vol. 100. pp. 635–640. 1978. 

–––––. “An Experimental Study and Theoretical Study of Natural Convection in the Annulus 

Between Horizontal Concentric Cylinders.” Journal of Fluid Mechanics. Vol. 74, Part 4. pp. 

695–719. 1976. 

Manepally, C., S. Green., F. Viana, and R. Fedors. “Evaluation of In-Drift Heat Transfer 

Processes.” San Antonio, Texas: CNWRA. 2006. 

Menter, F.R. “Two-Equation Eddy-Viscosity Turbulence Models for Engineering Applications.” 

American Institute of Aeronautics and Astronautics Journal. Vol. 32, No. 8. pp. 1,598–1,605. 

1994. 

Popiel, C.O. and J. Wojtkowiak. “Simple Formulas for Thermophysical Properties of Liquid 

Water for Heat Transfer Calculations (0oC to 150oC).” Heat Transfer Engineering. Vol. 19, 

No. 3. pp. 87–101. 1998. 

Reynolds, W.C. “Thermodynamic Properties in SI: Graphs, Tables, and Computational 

Equations for Forty Substances.” Palo Alto, California: Stanford University, Department of 

Mechanical Engineering. 1979. 

Rhie C.M. and W.L. Chow. “Numerical Study of the Turbulent Flow Past an Airfoil With Trailing 

Edge Separation.” American Institute of Aeronautics and Astronautics Journal. Vol. 21. 

pp. 1,523–1,532. 1983. 

Salari, K. “Adaptive Coupling Techniques for CFD and Porous Flow.” Engineering Research and 

Technology Report. Livermore, California: Lawrence Livermore National Laboratory. 2005 

Sparrow, E.M., W.J. Minkowycz, and M. Saddy. “Forced Convection Condensation in the 

Presence of Noncondensables and Interfacial Resistance.” International Journal—Heat Mass 

Transfer. Vol. 10. pp. 1,829–1,845. 1967. 

Van Doormal, J.P. and G.D. Raithby. “Enhancements of the SIMPLE Method for Predicting 

Incompressible Fluid Flows.” Numerical Heat Transfer. Vol. 7. pp. 147–163. 1984.