Thermal Analysis of 

Horizontal Storage 

Casks for Extended 

Storage Applications

Office of Nuclear Regulatory Research

NUREG/CR-7191
NRC Reference Material 

As of November 1999, you may electronically access 

NUREG-series publications and other NRC records at 

NRC’s Public Electronic Reading Room at 

http://www.nrc.gov/reading-rm.html. Publicly released 

records include, to name a few, NUREG-series 

publications; Federal Register notices; applicant, 

licensee, and vendor documents and correspondence; 

NRC correspondence and internal memoranda; bulletins 

and information notices; inspection and investigative 

reports; licensee event reports; and Commission papers 

and their attachments. 

NRC publications in the NUREG series, NRC 

regulations, and Title 10, “Energy,” in the Code of

Federal Regulations may also be purchased from one 

of these two sources. 

1. The Superintendent of Documents

 U.S. Government Printing Office

Mail Stop SSOP

Washington, DC 20402–0001

Internet: bookstore.gpo.gov

Telephone: 202-512-1800

Fax: 202-512-2250 

2. The National Technical Information Service 

Springfield, VA 22161–0002

www.ntis.gov

1–800–553–6847 or, locally, 703–605–6000 

A single copy of each NRC draft report for comment is 

available free, to the extent of supply, upon written 

request as follows:

Address: U.S. Nuclear Regulatory Commission

 Office of Administration 

Publications Branch

 Washington, DC 20555-0001

E-mail: DISTRIBUTION.RESOURCE@NRC.GOV

Facsimile: 301–415–2289 

Some publications in the NUREG series that are 

posted at NRC’s Web site address 

http://www.nrc.gov/reading-rm/doc-collections/nuregs 

are updated periodically and may differ from the last 

printed version. Although references to material found on 

a Web site bear the date the material was accessed, the 

material available on the date cited may subsequently be 

removed from the site.

Non-NRC Reference Material 

Documents available from public and special technical 

libraries include all open literature items, such as books, 

journal articles, transactions, Federal Register notices, 

Federal and State legislation, and congressional reports. 

Such documents as theses, dissertations, foreign reports 

and translations, and non-NRC conference proceedings 

may be purchased from their sponsoring organization. 

Copies of industry codes and standards used in a 

substantive manner in the NRC regulatory process are 

maintained at—

The NRC Technical Library

Two White Flint North

11545 Rockville Pike

Rockville, MD 20852–2738 

These standards are available in the library for reference 

use by the public. Codes and standards are usually 

copyrighted and may be purchased from the originating 

organization or, if they are American National Standards, 

from—

American National Standards Institute

11 West 42nd Street

New York, NY 10036–8002

www.ansi.org

212–642–4900 

AVAILABILITY OF REFERENCE MATERIALS

IN NRC PUBLICATIONS

Legally binding regulatory requirements are stated only 

in laws; NRC regulations; licenses, including technical 

specifications; or orders, not in NUREG-series 

publications. The views expressed in contractorprepared publications in this series are not necessarily 

those of the NRC.

The NUREG series comprises (1) technical and 

administrative reports and books prepared by the staff 

(NUREG–XXXX) or agency contractors (NUREG/CR–

XXXX), (2) proceedings of conferences (NUREG/CP–

XXXX), (3) reports resulting from international 

agreements (NUREG/IA–XXXX), (4) brochures 

(NUREG/BR–XXXX), and (5) compilations of legal 

decisions and orders of the Commission and Atomic and 

Safety Licensing Boards and of Directors’ decisions

under Section 2.206 of NRC’s regulations (NUREG–

0750).

DISCLAIMER: This report was prepared as an account 

of work sponsored by an agency of the U.S. 

Government. Neither the U.S. Government nor any 

agency thereof, nor any employee, makes any warranty, 

expressed or implied, or assumes any legal liability or 

responsibility for any third party’s use, or the results of 

such use, of any information, apparatus, product, or 

process disclosed in this publication, or represents that 

its use by such third party would not infringe privately 

owned rights.
Thermal Analysis of 

Horizontal Storage 

Casks for Extended 

Storage Applications

Manuscript Completed: February 2014

Date Published: December 2014

Prepared by: 

Kaushik Das, Debashis Basu, and Gary Walter

Center for Nuclear Waste Regulatory Analyses 

Southwest Research Institute® 

6220 Culebra Road 

San Antonio, TX 78238-5166

S. Gonzalez, NRC Project Manager

A. Velazquez-Lozada, NRC Project Manager

G. Zigh, NRC Technical Lead

Office of Nuclear Regulatory Research

NUREG/CR-7191

 

iii

ABSTRACT

A computational fluid dynamics (CFD) model that can be used to reliably predict temperature 

distributions for long-term storage was developed for the ventilated horizontal dry storage

system containing pressurized water reactor (PWR) fuel at the Calvert Cliffs Independent Spent 

Fuel Storage Installation (ISFSI). The model was developed systematically by (i) constructing a 

three dimensional (3-D) model for two horizontal modules at this ISFSI, (ii) comparing model 

results with measured temperature data, (iii) performing a parametric analysis to assess the 

impact of model inputs on temperature predictions, (iv) using the developed model to predict 

temperature distributions of critical components for a storage period of 300 years, and 

(v) conducting numerical uncertainty analyses of the solutions using the grid convergence index 

(GCI) method. 

Three-dimensional (3-D) numerical models were developed using the commercial CFD software 

FLUENT, Version 14.5 [ANSYS Inc. 2012]. The horizontal storage modules (HSM) modeled in 

the analysis are the HSM–1 and HSM–15 units at Calvert Cliffs. Each of the storage modules 

contains a horizontal dry shielded canister (DSC) loaded with 24 fuel assemblies. The concrete 

storage units and DSC internal structures such as spacer plates, tie rods, fuel assemblies, and 

supporting rails were included in the model. Individual fuel rods were not explicitly modeled; 

instead, the volume within each fuel assembly was represented as a porous medium with 

specified frictional and inertial flow resistances to helium movement. An orthotropic 

temperature-dependent equivalent thermal conductivity was used to model conductive and 

radiative heat transfer in the porous zone. The discrete ordinate (DO) model was used for 

radiation, while the low Reynolds number k-ε model with buoyancy was used to 

simulate turbulence.

The baseline simulations were compared to measured temperature data for the HSM–1 and 

HSM–15 units from June 27 and 28, 2012 [Suffield et al. 2012]. The steady state temperature 

distributions on different components of the storage modules, and the corresponding air and 

helium circulation patterns, were calculated using ambient conditions representative of the 

measurement period. The models predicted a high-temperature region on the upper half of the 

curved surface of the DSC shell, resulting from a combination of temperature-driven natural 

convection flow of the helium coolant inside the DSC and buoyancy-driven upward flow of 

heated air around the DSC. Computed temperatures at selected locations on the DSC shell 

were compared with measured data. The simulated temperatures for both modules were, 

higher than the measured values, but the agreement between the modeled and measured data 

was better near the end cap of the DSC shell. Aspects of the temperature data collection 

suggest that there were a number of uncertainties involved with the method used for 

temperature measurement and the observed values were not representative of a normally 

operated closed module.

Sensitivity analyses of the simulated results included evaluations of the near wall mesh 

refinement, turbulence model, porous media resistance, insolation, and evaluation of the heat 

transfer coefficient using an extended domain analysis. The study showed that near wall mesh 

refinement above the baseline mesh had minimal impact on the temperature distribution.

Therefore, the baseline mesh was deemed adequate to perform the benchmark studies, other 

parametric analyses, and the 300-year thermal evolution study. The turbulence models

evaluated included the standard k-ε model, low Reynolds number k-ε model, standard SST k-ω

model, and SST k-ω model with buoyance effects. Among the studied turbulence models, 

simulated results with the low Reynolds number k-ε model predicted the lowest temperatures. 

The study of porous media resistance showed that as long as some flow resistance exists within 

the porous zone of the fuel assemblies, the temperature distribution is not very sensitive to the 
iv

specified porous media resistance; however, the fundamental pattern of temperature distribution 

is altered if resistance is completely eliminated from this zone. Accordingly, reasonable 

resistance coefficients are necessary in this zone in order to accurately capture the temperature 

field. The sensitivity study also showed that insolation impacts the temperatures calculated on 

the internal components by a few degrees. To obtain an independent verification of the heat 

transfer coefficient on the exposed concrete surface of the module, a simulation was performed 

using an extended domain that included the surrounding atmosphere. Results indicated that the 

heat transfer coefficients used for the baseline simulation were reasonable. 

The 300-year simulations for both the HSM–1 and HSM–15 configurations showed that (i) the 

drop in maximum cladding temperature along the profile is relatively rapid in the first 100 years 

and becomes more gradual thereafter and (ii) the minimum cladding temperature varies by a 

small amount as time progresses. All the component temperatures showed an asymptotic 

decline towards the ambient temperature over time.

The GCI method developed by Roache [1998, 2003] and adopted in the ASME standard V&V

20-2009 [Celik et al. 2008; ASME 2009] was used to calculate the observed order of accuracy 

and associated numerical uncertainty for the simulated results. Four levels of computational 

grids were used in the GCI analysis. The simulated results obtained from the different grid 

levels showed little variation. The apparent order of accuracy calculated using different 

combinations of solutions was not constant, which can be attributed to minor numerical 

oscillations in the results and possibly round-off errors. An alternative estimate of numerical 

uncertainty was made using an order of accuracy of 1 and a higher factor of safety. From this 

estimation, it was found that the GCI is less than 4 percent, indicating that there was no 

significant deviation of the predicted variables due to mesh refinement. 
v

TABLE OF CONTENTS

Section Page

ABSTRACT ........................................................................................................................... iii

TABLE OF CONTENTS..............................................................................................................v

LIST OF FIGURES ....................................................................................................................ix

LIST OF TABLES......................................................................................................................xv

EXECUTIVE SUMMARY .........................................................................................................xvii

ACKNOWLEDGMENTS........................................................................................................... xxi

ABREVIATIONS AND ACRONYMS........................................................................................xxiii

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

1.1 Regulatory and Technical Background..............................................................1-1

1.2 Brief Literature Survey ......................................................................................1-1

1.3 Objective and Scope of Work............................................................................1-3

1.4 Description of the Dry Cask Storage Modules (HSM–1 and HSM–15) .............1-3

2. NUMERICAL MODELING METHODS..........................................................................2-1

2.1 ANSYS-FLUENT Fluid Flow Solver...................................................................2-1

2.2 Governing Equations ........................................................................................2-1

2.2.1 Continuity and Momentum Equations ....................................................2-1

2.2.2 Energy Equation....................................................................................2-2

2.3 Solid Modeling and Mesh Generation for HSM–1 and HSM–15 .......................2-3

2.3.1 Grid Generation.....................................................................................2-4

2.4 Boundary and Operating Conditions .................................................................2-7

2.5 Numerics, Radiation Heat Transfer, and Turbulence Models for 

the Simulations ...............................................................................................2-12

2.6 Fuel Assembly Model......................................................................................2-13

2.6.1 Porous Media Calculation....................................................................2-13

2.6.2 Effective Thermal Conductivity Calculation..........................................2-13

2.6.3 Fuel Assembly Loading .......................................................................2-14

2.7 Material Properties Used in the Simulations....................................................2-14

3. BENCHMARK MODELING AND SENSITIVITY ANALYSES........................................3-1

3.1 Description of Field Temperature Measurement................................................3-1

3.2 Results of Baseline Simulations ........................................................................3-1

3.2.1 HSM–1 ..................................................................................................3-1

3.2.1.1 Temperature and Velocity Patterns.....................................3-2

3.2.1.2 Temperature Line Probes ...................................................3-4

3.2.1.3 Comparison With Measured Data .....................................3-11

3.2.2 HSM–15 ..............................................................................................3-14

3.2.2.1 Temperature Field.............................................................3-14

3.2.2.2 Velocity Field ....................................................................3-19

3.2.2.3 Comparison with Measured Data......................................3-19

3.3 Sensitivity Studies...........................................................................................3-19

3.3.1 Effect of Mesh Refinement ..................................................................3-21

3.3.2 Effect of Turbulence Models................................................................3-22

3.3.3 Effect of Porous Media Resistance......................................................3-27

3.3.4 Effect of Insolation...............................................................................3-31

3.3.5 Effect of Including the Surrounding Atmosphere in the 

Computational Domain........................................................................3-33
vi

TABLE OF CONTENTS

(continued)

Section Page

3.4 Observations From Benchmark Study.............................................................3-36

4. THERMAL EVOLUTION SIMULATIONS......................................................................4-1

4.1 Numerical Approach to Thermal Evolution Analysis ..........................................4-1

4.2 Results of Thermal Evolution Simulations .........................................................4-1

4.3 Observations from the Thermal Evolution Study .............................................4-16

5. UNCERTAINTY QUANTIFICATION AND GRID CONVERGENCE INDEX (GCI) 

ANALYSIS....................................................................................................................5-1

5.1 Uncertainty Analysis: Definitions and Concepts ...............................................5-1

5.1.1 Code Validation.....................................................................................5-1

5.1.2 Code Verification...................................................................................5-2

5.1.3 Solution Verification...............................................................................5-2

5.2 Grid Convergence Index Calculation Method ....................................................5-3

5.3 Uncertainty Quantification and GCI Analysis for HSM–1 Module ......................5-5

5.3.1 Assessment of Asymptotic Behavior of the Solution ..............................5-9

5.3.2 Iteration Convergence .........................................................................5-11

5.3.3 Other Numerical Issues.......................................................................5-11

5.3.4 Modified Estimation of GCI for HSM-1.................................................5-14

5.4 Uncertainty Quantification and GCI Analysis for HSM–15 Module ..................5-14

5.5 Observations from GCI Analysis .....................................................................5-18

6. SUMMARY AND CONCLUSIONS................................................................................6-1

6.1 Baseline Simulations and Comparison of Simulated to 

Measured Temperatures...................................................................................6-1

6.2 Sensitivity Analyses ..........................................................................................6-2

6.3 Simulation of Long-Term Storage Conditions....................................................6-3

6.4 Uncertainty Analysis .........................................................................................6-3

7. REFERENCES.............................................................................................................7-1

APPENDIX A: NATURAL CONVECTION CORRELATION AND MODEL FOR VENT

SCREENS AT INLET..................................................................................... A–1

A.1 Natural Convection Correlation....................................................................... A–1

A.1.1 Correlations for Vertical Flat Plates ..................................................... A–2

A.1.2 Correlation for Horizontal Flat Plates................................................... A–3

A.2 Modeling of Vents Screens at Inlet.................................................................. A–3

APPENDIX B: MATERIAL PROPERTIES USED IN THE SIMULATIONS.............................. B–1

APPENDIX C: POROUS MEDIA MODEL .............................................................................. C–1
vii

TABLE OF CONTENTS

(continued)

Section Page

APPENDIX D: CALCULATION OF FUEL EFFECTIVE THERMAL CONDUCTIVITY ............. D–1

D.1 Radial Thermal Conductivity ........................................................................... D–1

D.2 Axial Thermal Conductivity.............................................................................. D–3

APPENDIX E: FUEL ASSEMBLY LOADING.......................................................................... E–1

ix

LIST OF FIGURES

Figure Page

1-1 NUHOMS horizontal storage modules in Calvert Cliffs Nuclear Power 

Station ISFSI ................................................................................................................1-4

1-2 The concrete storage module and DSC shell with support rail for unit HSM–1 ............1-5

1-3 DSC sell and internal components: spacer plates, fuel assemblies and tie rods 

for units HSM–1............................................................................................................1-5

2-1 Aerial view of the Calvert Cliffs Nuclear Power Station ISFSI housing HSM–1

and HSM–15 modules ..................................................................................................2-3

2-2 Schematics showing the geometry of the DSC.............................................................2-4

2-3 Three-dimensional schematic view of the HSM–1 module............................................2-5

2-4 Three-dimensional view of the DSC internal showing the support rod, guide

sleeve, spacer disk, and guide rails ..............................................................................2-5

2-5 Three-dimensional schematic for the fuel assemblies, tie rods, and spacer

plates within the DSC shell...........................................................................................2-6

2-6 Two-dimensional cross-sectional view showing the computational grid

for the fuel region .........................................................................................................2-6

2-7 Three-dimensional view of the computational grid (exterior view) of the

HSM–1 module.............................................................................................................2-8

2-8 Transverse slice through center line showing computational grid of HSM–1,

highlighting the airflow path ..........................................................................................2-8

2-9 Planar slice through mid-section of the DSC showing computational mesh

of HSM–1 module in transverse direction .....................................................................2-9

2-10 Enlarged view of the computational grid in the guide sleeve, support rod, 

and spacer plate region ................................................................................................2-9

2-11 Three-dimensional view of the computational grid for the DSC internal and

guide rails...................................................................................................................2-10

2-12 Three-dimensional view of the DSC internal showing the fuel region,

support rod, guide sleeve, spacer disk, and guide rails...............................................2-10

2-13 Fuel assembly loading identification structure for HSM–1 module ..............................2-15

2-14 Fuel assembly loading identification structure for HSM–15 module ............................2-15

3-1 Axial and transverse cross-sectional planes used to illustrate temperature

and velocity distributions ..............................................................................................3-2

3-2 Contours of DSC shell surface temperature for HSM–1 baseline case .........................3-3

3-3 Temperature distributions for the fuel region, spacer plates, and support rods for

HSM–1 baseline case...................................................................................................3-3

3-4 Air and helium temperature distributions along transverse cross section

for HSM–1 baseline case..............................................................................................3-5

3-5 Temperature distribution for the solid region at transverse cross section

for HSM–1 baseline case..............................................................................................3-5

3-6 Temperature distribution for the fluid region along the axial cross-sectional plane

for HSM–1 baseline case..............................................................................................3-6

3-7 Temperature distribution for the solid region along the axial cross-sectional plane

the HSM–1 module.......................................................................................................3-6

3-8 Streamlines showing the airflow path and airflow circulation within the 

HSM–1 module.............................................................................................................3-7

3-9 Velocity vectors at the mid-plane in the transverse cross section for HSM–1................3-7
x

LIST OF FIGURES

(continued)

Figure Page

3-10 Velocity vectors at the axial cross section for HSM–1 module ......................................3-8

3-11 Line probe locations for extracting temperature distribution along the axial 

direction: locations shown on the DSC shell and on the heat shield ............................3-8

3-12 Line probe locations for extracting temperature distribution along the axial

direction in the fuel region: locations shown for fuel assemblies 1 and 23 ...................3-9

3-13 Computed temperature distributions along the axial distance for top of DSC

shell and right of DSC shell ........................................................................................3-10

3-14 Computed temperature distributions along the axial distance for 

top of heat shield and right of heat shield ...................................................................3-10

3-15 Computed temperature distributions along the axial distance for fuel region: 

distribution shown for fuel assemblies 1 and 23 for HSM–1 baseline case .................3-11

3-16 Thermocouple locations shown at 0.0-, 0.51-, and 1.02-m axial distance ...................3-12

3-17 Thermocouple locations shown at the top, right side, left side, rail right, and

rail left .......................................................................................................................3-13

3-18 Contours of DSC shell surface temperature for HSM–15 baseline case .....................3-15

3-19 Temperature distributions for the fuel assembly, spacer plates, and 

support rod for HSM–15 baseline case.......................................................................3-15

3-20 Air and helium temperature distributions along transverse cross section 

for HSM–15 baseline case..........................................................................................3-16

3-21 Solid region (concrete, spacer plates) temperature distributions at the 

transverse cross section for HSM–15 baseline case....................................................316

3-22 Temperature distribution for air and helium at axial cross-sectional plane

for HSM–15 baseline case..........................................................................................3-17

3-23 Temperature distribution for the solid regions at axial cross-sectional plane

for HSM–15 baseline case..........................................................................................3-17

3-24 Computed temperature distributions along the axial distance for top and

right of DSC shell for HSM–15 module .......................................................................3-18

3-25 Computed temperature distributions along the axial distance for top of heat

shield and right of heat shield .....................................................................................3-18

3-26 Computed temperature distributions along the axial distance for fuel region:

distribution shown for fuel assemblies 1 and 23 for HSM–15 baseline case ...............3-19

3-27 Streamlines showing the airflow path and airflow circulation within the 

HSM–15 module.........................................................................................................3-20

3-28 Velocity vector distribution at the transverse cross section for HSM–15 module.........3-20

3-29 Velocity distribution at the axial cross section for HSM–15 module.............................3-21

3-30 Computed temperature distributions along the axial distance for right side of

DSC shell: mesh refinement for HSM–1 module........................................................3-23

3-31 Computed temperature distributions along the axial distance for top of 

DSC shell: mesh refinement for HSM–1 module........................................................3-23

3-32 Computed temperature distributions along the axial distance for fuel assembly 1:

mesh refinement for HSM–1 module ..........................................................................3-24

3-33 Computed temperature distributions along the axial distance for right side of 

DSC shell: effect of turbulence model for HSM–1 module .........................................3-25

3-34 Computed temperature distributions along the axial distance for top of

DSC shell: effect of turbulence models for HSM–1 module........................................3-26
xi

LIST OF FIGURES

(continued)

Figure Page

3-35 Computed temperature distributions along the axial distance for fuel assembly 1: 

effect of turbulence model for HSM–1 module............................................................3-26

3-36 Computed temperature distributions along the axial distance for fuel 

assembly 23: effect of turbulence model for HSM–1 module .....................................3-27

3-37 Computed temperature distributions along the axial distance for right side of 

DSC shell: effect of porous media for HSM–1 module ...............................................3-28

3-38 Computed temperature distributions along the axial distance for top of DSC 

shell: effect of porous media for HSM–1 module .......................................................3-29

3-39 Computed temperature distributions along the axial distance for fuel assembly 1:

 effect of porous media for HSM–1 module..................................................................3-29

3-40 Computed temperature distributions along the axial distance for fuel 

assembly 23: effect of porous media for HSM–1 module...........................................3-30

3-41 Computed temperature distributions along the axial distance for right side and

top of DSC shell: effect of insolation for HSM–1 module............................................3-32

3-42 Computed temperature distributions along the axial distance for fuel 

assemblies 1 and 23: effect of insolation for HSM–1 module.....................................3-32

3-43 Extended domain for the study of the effect of atmosphere for HSM–1 module..........3-34

3-44 Temperature contours at an axial mid-plane of the extended domain: effect of 

atmosphere for HSM–1 module..................................................................................3-35

3-45 Velocity contours at an axial mid-plane of the extended domain: effect of 

atmosphere for HSM–1 module..................................................................................3-35

3-46 Computed heat transfer coefficient distribution over the HSM–1 module: 

effect of extended domain and atmosphere................................................................3-36

4-1 Variation of maximum and minimum cladding temperature over 300 years for

HSM–1 module.............................................................................................................4-5

4-2 Variation of maximum and minimum fuel sleeve temperature over 300 years for

HSM–1 module.............................................................................................................4-5

4-3 Variation of maximum and minimum spacer plates temperature over 300 years

for HSM–1 module........................................................................................................4-6

4-4 Variation of maximum and minimum DSC shell temperature over 300 years for

HSM–1 module.............................................................................................................4-6

4-5 Variation of maximum and minimum heat shield temperature over 300 years for

HSM–1 module.............................................................................................................4-7

4-6 Air and helium temperature distributions at a transverse plane for 

HSM–1: year 1993.......................................................................................................4-7

4-7 Air and helium temperature distributions at a transverse plane for 

HSM–1: year 2093 (100 years after storage)...............................................................4-8

4-8 Air and helium temperature distributions at a transverse plane for 

HSM–1: year 2193 (200 years after storage)...............................................................4-8

4-9 Air and helium temperature distributions at a transverse plane for 

HSM–1: year 2193 (300 years after storage)...............................................................4-9
xii

LIST OF FIGURES

(continued)

Figure Page

4-10 Variation of maximum and minimum cladding temperature over 300 years 

for HSM–15 module......................................................................................................4-9

4-11 Variation of maximum and minimum fuel sleeve temperature over 300 years

for HSM–15 module....................................................................................................4-10

4-12 Variation of maximum and minimum spacer plates temperature over 300 years

for HSM–15 module....................................................................................................4-10

4-13 Variation of maximum and minimum DSC shell temperature over 301 years 

for HSM–15 module....................................................................................................4-11

4-14 Variation of maximum and minimum heat shield temperature over 301 years

for HSM–15 module....................................................................................................4-11

4-15 Air and helium temperature distributions at a transverse plane for 

HSM–15: year 1996...................................................................................................4-12

4-16 Air and helium temperature distributions at a transverse plane for 

HSM–15: year 2097 (101 years after storage)...........................................................4-12

4-17 Air and helium temperature distributions at a transverse plane for 

HSM–15: year 2197 (201 years after storage)...........................................................4-13

4-18 Air and helium temperature distributions at a transverse plane for 

HSM–15: year 2297 (301 years after storage)...........................................................4-13

5-1 Temperature contours along the mid-plane for four mesh levels for 

simulations of HSM–1 at year 1993 ..............................................................................5-7

5-2 Deviation of temperature values calculated at different mesh levels for

HSM–1: deviations are calculated as percentage of fine mesh temperatures ..............5-8

5-3 Calculated apparent order using different mesh levels for HSM–1..............................5-10

5-4 Distribution of discretization error with mesh refinement for HSM–1 ...........................5-11

5-5 Variation of convergence coefficient with mesh refinement for HSM–1.......................5-12

5-6 Variation of average temperature of selected monitor points with iteration,

showing overall convergence .....................................................................................5-12

5-7 Variation of average temperature of selected monitor points with iteration

for last 6,000 iterations, highlighting minor oscillations ...............................................5-13

5-8 Deviation of temperature values calculated at different mesh levels for

HSM–15: deviations are calculated as percentage of fine mesh temperatures ..........5-16

5-9 Calculated apparent order using different mesh levels for HSM–15............................5-17

5-10 Distribution of discretization error with mesh refinement for HSM–15 .........................5-17

5-11 Variation of convergence coefficient with mesh refinement for HSM–15.....................5-19

C–1 Longitudinal full view of the computational domain for the CE 14 × 14 fuel 

assembly for porous media parameter calculation showing the fuel rod ..................... C–2

C–2 1/4th cross-sectional domain for the CE 14 × 14 fuel assembly for porous

media parameter calculation....................................................................................... C–2

C–3 Two-dimensional view of the computational domain for the CE 14 × 14 fuel

assembly for porous media parameter calculation...................................................... C–3

D–1 Computational domain and grid for fuel, clad, and helium gas grids in the

effective thermal conductivity calculation .................................................................... D–2

D–2 Predicted temperature contours for an edge temperature of 300 K [80.3 °F] .............. D–2
xiii

LIST OF FIGURES

(continued)

Figure Page

E–1 Fuel assembly loading identification structure for HSM–1 module .............................. E–2

E–2 Fuel assembly loading identification structure for HSM–15 module ............................ E–2

E–3 Axial distribution of the peaking factor used in current simulations ............................. E–4

xv

LIST OF TABLES

Table Page

2-1 Total number of cells for different grid levels used in the GCI study............................2-11

3-1 Comparison of measured and predicted temperatures for the baseline 

case simulations for HSM–1 module ..........................................................................3-12

3-2 Computational mesh used in the mesh refinement study............................................3-24

3-3 Comparison of predicted temperatures for the three different levels of 

mesh: simulations for HSM–1 module .......................................................................3-24

3-4 Comparison of predicted temperatures for different turbulence models:

simulations for HSM–1 module...................................................................................3-28

3-5 Resistance coefficients used porous media resistance parametric study....................3-30

3-6 Comparison of predicted temperatures (K) for different porous media resistance

values: simulations for HSM–1 module......................................................................3-31

3-7 Comparison of predicted temperatures (K) for with and without

insolation: simulations for HSM–1 module .................................................................3-34

4-1 Calculated decay heat load values at different year for HSM–1 module .......................4-2

4-2 Calculated decay heat load values at different year for HSM–15 module .....................4-3

4-3 Calculated temperature values at different years for HSM–1 module .........................4-14

4-4 Calculated temperature values at different years for HSM–15 module........................4-15

5-1 Mesh levels used in GCI analysis of HSM–1 module....................................................5-6

5-2 Computed maximum and minimum cladding temperature for HSM–1

module (1993) ..............................................................................................................5-8

5-3 Computed maximum and minimum cladding temperature for HSM–1

module (2012) ..............................................................................................................5-8

5-4 Computed maximum and minimum cladding temperature for HSM–1

module (2093) ..............................................................................................................5-8

5-5 Apparent order (p) and calculated GCI for HSM–1 module...........................................5-9

5-6 GCI value for a conservative bounded value of p = 1 for HSM–1 module ...................5-15

5-7 Mesh levels used in GCI analysis of HSM–15 module................................................5-15

5-8 Computed maximum and minimum cladding temperature for HSM–15 

module (1996) ............................................................................................................5-15

5-9 Computed maximum and minimum cladding temperature for HSM–15 

module (2012) ............................................................................................................5-15

5-10 Computed maximum and minimum cladding temperature for HSM–15 

module (2097) ............................................................................................................5-15

5-11 Apparent order (p) and calculated GCI for HSM–15 module.......................................5-17

5-12 GCI value for a conservative bounded value of p = 1 for HSM–15 module .................5-19

B–1 Properties of air used in the simulations ..................................................................... B–1

B–2 Properties of helium used in the simulations............................................................... B–1

B–3 Properties of concrete used in the simulations............................................................ B–1

B–4 Properties of carbon steel used in the simulations...................................................... B–1

B–5 Properties of stainless steel used in the simulations................................................... B–2

B–6 Properties of lead used in the simulations .................................................................. B–2

B–7 Emissivity of materials used in the simulations ........................................................... B–2
xvi

LIST OF TABLES

(continued)

Table Page

C–1 Dimensions used in the simulations............................................................................ C–3

C–2 Inlet condition and parameters for calculations of loss factors .................................... C–4

D–1 Keff of fuel region in the radial direction with helium backfill......................................... D–3

D–2 Keff of fuel region in the axial direction with helium backfill .......................................... D–3

E–1 Calculated decay heat load values at year 2012 for HSM–1 module .......................... E–3

E–2 Calculated decay heat load values at year 2012 for HSM–15 module ........................ E–3
xvii

EXECUTIVE SUMMARY

The U.S. Nuclear Regulatory Commission (NRC) regulates storage of spent nuclear fuel (SNF) 

from commercial nuclear power plants. An increasing amount of the SNF is in dry storage 

systems, mostly at current and decommissioned plants. The existing regulatory framework in 

the Code of Federal Regulations (CFR) for storage (10 CFR 72) and transportation (10 CFR 71)

of SNF supports dry cask storage for at least 80 years (i.e., a 40-year initial licensing term, 

followed by a license renewal for a term of up to 40 years, although many of the existing 

facilities were licensed for an initial term of 20 years under the regulations in place at the time). 

The NRC staff is currently considering the performance of storage systems over an initial

300-year period following discharge of the SNF from the reactor, in advance of potential 

changes to the existing regulatory framework. 

This report describes an independent computational analysis designed to predict the thermal 

behavior of dry cask storage systems for periods of up to 300 years. The simulations were 

performed for ventilated dry cask horizontal storage modules HSM–1 and HSM–15 at the 

Calvert Cliffs Independent Spent Fuel Storage Installation (ISFSI). Three-dimensional (3-D) 

numerical simulations were performed using the computational fluid dynamics (CFD) software 

ANSYS-FLUENT® Version 14.5 [ANSYS 2012]. Initially, a baseline model of the storage 

modules was developed, and the results highlighted the general pattern of the flow and 

temperature fields within the storage units. Temperature values calculated from the baseline 

study were compared with available measured data. Subsequently, a sensitivity study was 

conducted to understand the effects on the numerical solutions of different modeling 

parameters, such as the turbulence model, porous media resistance, and insolation. The 

baseline computational model was then used to predict the thermal behavior of the dry storage 

modules over 300 years of storage. Finally, a study to quantify numerical uncertainty was 

performed using the grid convergence index (GCI) method.

The HSM–1 and HSM–15 ventilated dry cask storage modules each contain a dry shield 

canister (DSC) that is loaded horizontally. The DSC shell holds 24 fuel assemblies, each 

containing an arrangement of Combustion Engineering (CE) 14 × 14 fuel rods. The void spaces 

within the DSC shell contain helium as a coolant, maintained above atmospheric pressure. 

Ambient air enters the base of each storage module, passes over the surface of the DSC shell, 

and exits from vents at the top of the module. The 3-D computational model developed in this 

study explicitly includes the concrete storage units and DSC shells, along with internal 

structures such as spacer plates, tie rods, fuel assemblies, and the supporting rails. The 

individual fuel rods within an assembly were not explicitly represented; instead, the volume 

within each fuel assembly was modeled as a porous medium to represent flow of the helium 

coolant. Inertial and frictional flow resistance factors, which were obtained using a separate 3-D 

CFD model of the fuel assembly, were specified for helium movement in the porous zone as 

shown in Appendix C. An orthotropic temperature-dependent equivalent thermal conductivity, 

which was calculated based on the CE 14 × 14 fuel configuration, was used to model the 

conductive and radiative heat transfer in the porous zone. Hexahedral grids were used to mesh 

the fluid and solid volumes in the domain. The discrete ordinate (DO) model was used for 

radiation. The low Reynolds number k-ε model with buoyancy effect consideration was used to 

simulate turbulence.

All of the simulations used ambient conditions corresponding to June 27 and 28, 2012, when 

HSM–1 and HSM–15 were opened for temperature measurements. The baseline simulations 

calculate the steady-state temperature distribution for different components of the storage 

module and the circulation pattern for ventilation air and the helium coolant. In general, a 
Executive Summary

xviii

high-temperature region resides on the upper half of the curved surface of the DSC shell, 

resulting from temperature-driven natural convection flow of the helium coolant inside the DSC

and the flow of ventilation air around the DSC. Computed temperature values at selected 

locations on the DSC shell surface were compared with measured data. The computed results 

were consistently higher than the measured temperatures. Measurement uncertainty may be 

significant because of the difficult access conditions and exchange of air through the open DSC 

door, and may have contributed to the differences between computed and measured values.

Once the 3-D baseline model was developed and benchmarked with experimental data, it was 

used to perform a series of sensitivity studies. The sensitivity studies involved analyzing the 

effects of near wall mesh refinement, turbulence model, porous media resistance, an extended 

computational domain, and insolation. The mesh refinement study was done by adapting and 

refining the mesh near the DSC shell wall. A comparison of temperature distributions among

solutions using three levels of adapted mesh showed minimal impacts of near wall refinement. 

Therefore, the baseline mesh was considered adequate to perform the present analyses. The 

turbulence models considered in the sensitivity analysis included the baseline low Reynolds 

number k-ε model, the standard k-ε model with scalable wall function and buoyancy effects, the 

shear stress transport (SST) k-ω model, and the SST k-ω model with buoyancy effects. A 

comparison of results obtained using the different turbulence models showed that the low

Reynolds number k-ε model predicted the lowest temperatures. To study the sensitivity to the 

specified porous media resistance, the resistances were increased and decreased by an order 

of magnitude from the baseline values. It was found that varying the resistances compared to 

baseline values did not affect the temperature distribution. However, case studies with zero 

resistance in the porous media and with no porous media assumption predicted a fundamentally 

altered temperature distribution pattern, indicating that a reasonable specification of porous 

resistance is necessary to obtain meaningful results. For the sensitivity analysis on the effect of 

insolation, a simulation was performed with no insolation flux at the boundary wall. Calculated 

temperatures on the DSC shell with no insolation were found to be lower than the baseline 

simulation results by a few degrees, indicating that insolation minimally impacts the internal 

component temperatures within the storage system. An analysis was performed to assess the 

validity of the convective heat transfer coefficient specified at the outer exposed boundary walls 

by extending the computational domain to include the surrounding atmosphere. It was found 

that the average convective heat transfer coefficient calculated from the extended domain 

calculation was comparable to the specified value at the boundaries of the baseline case.

The baseline computational model was used to predict the long-term (300 years) thermal 

behavior of the dry storage casks. Multiple steady-state cases were simulated at discrete time 

points. For the first 100 years, steady-state simulations were performed at 5-year intervals, 

dropping to 10-year intervals thereafter. In both the HSM–1 and HSM–15 configurations, (i) the 

drop in maximum cladding temperature along the profile is relatively rapid in the first 100 years 

and becomes more gradual thereafter and (ii) the minimum cladding temperature drops to a few 

degrees above the ambient temperature within 100 years. All of the component temperatures 

show an asymptotic decrease with time.

The GCI method, applied to evaluate numerical uncertainty in the simulations, uses the 

Richardson extrapolation approach to quantify the discretization error of the computed solution. 

Four levels of computational grids were used in the analysis: the baseline grid, a coarser grid, 

and two levels of finer grids. Maximum and minimum cladding temperatures were used as

target variables to estimate the uncertainty of the numerical simulations. These variables were 

selected as representative parameters because other component temperatures follow the same 

pattern over time. The calculated results obtained from different grid levels showed almost no 

variation. The minor variation between values could be attributed to effects such as machine 
Executive Summary

xix

round-off error. The apparent order of accuracy calculated using different combinations of 

solutions was not constant, and was larger than the theoretical order of accuracy in many of the 

cases. It appears that the small differences in temperature among different mesh levels 

disproportionally affected some terms associated with the calculation of apparent order, yielding

an artificially high value. The calculation of order p was also affected by minor oscillations of the 

results, even when overall convergence was achieved. To provide an estimate of numerical 

uncertainty, a theoretical order of accuracy of p = 1 was used with a higher factor of safety of 3. 

It was found that the GCI is less than 4 percent, indicating that there was no significant deviation 

due to mesh refinement.

xxi

ACKNOWLEDGMENTS

The authors gratefully acknowledge S. Stothoff for his technical review, G. Wittmeyer and 

D. Pickett for programmatic reviews, L. Mulverhill for her editorial review, A. Ramos for his 

administrative support, and G. Zigh and J. Solis for technical guidance. ANSYS Inc. is also 

acknowledged for its advice and help in ANSYS-FLUENT related issues used in this study. 

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. 

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

ANSYS® FLUENTTM Version 14.5 [ANSYS Inc. 2012a] was used to generate results for 

this report and is controlled in accordance with the CNWRA Technical Operating Procedure 

(TOP)–018, Development and Control of Scientific and Engineering Software. Post-processing 

was done using TECPLOT® 360TM Version 2009 [TECPLOT Inc. 2009] and ANSYS CFDPOST

[ANSYS 2012b].

xxiii

ABBREVIATIONS AND ACRONYMS

2-D two-dimensional

3-D three-dimensional

ASME American Society of Mechanical Engineers

CE Combustion Engineering

CFD computational fluid dynamics

CFR Code of Federal Regulations

CNWRA Center for Nuclear Waste Regulatory Analyses

DO discrete ordinate

DOE U.S. Department of Energy

DSC dry storage cask

EPRI Electric Power Research Institute

GCI grid convergence index

HSM horizontal storage module

ISFSI Independent Spent Fuel Storage Installation

LES large eddy simulation

MMS method of manufactured solution

NMSS Office of Nuclear Material Safety And Safeguards

NRC U.S. Nuclear Regulatory Commission

RES Office of Nuclear Regulatory Research

ORNL Oak Ridge National Laboratory

PDE partial differential equations

PNNL Pacific Northwest National Laboratory

PWR pressurized water reactor

RE Richardson extrapolation

RMS root mean square

RTD resistance temperature detector

RTE radiative transfer equation

SFS spent fuel storage

SNF spent nuclear fuel

SSC structures, systems, and components

SST shear stress transport

V&V verification and validation

1-1

1. INTRODUCTION

1.1 Regulatory and Technical Background

The U.S. Nuclear Regulatory Commission (NRC) regulates storage of spent nuclear fuel (SNF) 

from commercial nuclear power plants. The existing NRC regulatory framework for storage and 

transportation of SNF is specified in Title 10 of the Code of Federal Regulations (CFR) Part 72, 

“Licensing Requirements for the Independent Storage of Spent Nuclear Fuel, High-Level 

Radioactive Waste, and Reactor-Related Greater than Class C Waste” and 10 CFR Part 71, 

“Packaging and Transportation of Radioactive Material,” respectively. The current regulatory 

framework supports at least the first 80 years of dry cask storage (i.e., a 40-year initial licensing 

term, followed by a license renewal for a term of up to 40 years), although many of the existing 

facilities were licensed for an initial term of 20 years under the regulations in place at the time. 

As directed by the Commission [in SRM–COMSECY–10–0007, December 6, 2010]

1 in 

expectation of continued use of dry storage for extended periods of time, the NRC staff is 

examining the technical needs and potential changes to the regulatory framework that may be 

needed to continue licensing SNF storage over periods beyond 120 years. For this evaluation, 

the NRC staff has considered performance of the systems over an initial 300-year period 

following discharge of the SNF from the reactor. Modeling of long-term thermal behavior 

patterns forms an integral part of this assessment. Improvements in thermal calculations can 

enhance understanding of both the interior (e.g., fuel and cladding) and exterior (susceptibility of 

canisters to stress corrosion cracking) performance of system components. There is, therefore,

a technical need to assess the thermal behavior of spent fuel canisters for extended periods

of storage.

As part of this dry storage evaluation, the NRC Office of Nuclear Material Safety and 

Safeguards (NMSS) requested support from the NRC Office of Nuclear Regulatory Research 

(RES) to address the issues related to determining axial and radial temperature distributions of 

the fuel cladding as a function of time. This task involves performing state-of-the-art modeling 

of the thermal behavior of fuel in dry cask storage systems for periods of up to 300 years

[NRC 2012a]. This research will support risk-informed long-term storage regulations. For the 

foreseeable future, the staff anticipates dry cask storage will be the industry-preferred 

alternative for providing (i) additional storage capacity at operating and new reactors and (ii) the 

spent fuel management needs at decommissioned sites or at away-from-reactor sites. 

The current work by the Center for Nuclear Waste Regulatory Analyses (CNWRA®) provides 

technical assistance to the Division of Systems Analysis in RES to support the NRC evaluation 

of the long-term dry storage of SNF in the area of thermal modeling and uncertainty estimation 

of thermal models. This activity will help the CNWRA and NRC staffs prepare for potential 

technical issues related to the thermal behavior of SNF in storage.

1.2 Brief Literature Survey

Heat transfer in the cask is a complicated process because of the geometry of the spent fuel 

assemblies and dry cask canister, and the combined convection induced by the thermal 

radiation process [Heng et al. 2002]. Other factors that contribute to the overall heat transfer 

 

1

Vietti-Cook, A.L, “Staff Requirements–COMSECY–10–0007—Project Plan for Regulatory Program Review to 

Support Extended Storage and Transportation of Spent Nuclear,” email communication (December 6, 2010) to 

R.W. Borchardt, U.S. Nuclear Regulatory Commission, 2010.
Introduction

1-2

include the heat generation by the spent fuel, the ambient thermal boundary conditions, the 

filling medium within the storage cask, and the vertical or horizontal orientation of the cask. 

Some of the earliest studies performed in this area of cask modeling [Nishimura et al. 1996; 

Shibazaki et al. 1998] showed that the convective heat transfer in a storage cask depends on 

the Grashof number, the Prandtl number, and the relevant geometry parameters. Some of the 

earliest investigations [Arya and Keyhani 1990; Cannan and Klein 1998] focused on natural 

convection within the spent fuel assemblies. Gomez et al. [2005] performed two-dimensional 

(2-D) heat transfer simulations to predict the maximum cladding temperature in a boiling water 

reactor fuel assembly in a horizontal support basket. Simulated results showed that the 

predicted maximum cladding temperature is more sensitive to the assumed value of the fuel 

cladding emissivities when nitrogen is the backfill gas than when helium is used. Heng et al. 

[2002] found that with an increase in Rayleigh number, the dominant heat transfer mode in the 

cask changes from conduction to convection. Heng et al. [2002] also determined that at 

Rayleigh number = 1.3 × 109

, turbulent convective heat transfer is so strong that the 

temperature gradient is very steep near the wall of the cask, and that local convection 

circulation near the individual guide sleeves and the canister shell is more important to heat 

transfer than large-scale circulation within the canister. Walavalkar and Schowalter [2004]

performed a computational fluid dynamics (CFD) analysis of the VSC–17 spent fuel dry storage 

system using the FLUENT code with the thermal radiation equations solved for a 90° section of 

the VSC–17 system.

Zigh and Solis [2008] analyzed the effect of different turbulence models and boundary 

conditions on modeling uncertainty for the VSC–17 system. The simulations using the low 

Reynolds number k-ε and the transitional k-ω turbulence models compared well with measured 

temperatures. The simulation results also identified that the operating density and the ambient 

pressure significantly influence the results. The porous media model used to simulate the flow 

inside the storage canister in the current report was based on the approach described by 

Zigh and Solis [2008]. In previous numerical studies of dry storage canisters at CNWRA,

Das et al. [2010] simulated the VSC–17 system and compared the results with available 

experimental data. The simulated peak cladding temperature was slightly higher than the 

experimental data. However, the simulated results showed that the match between the 

computed results and the experimental data was better when nitrogen was used as a backfill 

gas inside the canister instead of helium.

Hoo et al. [2010] performed FLUENT numerical simulations of thermal hydraulic phenomena in 

a TN24P cask. They performed a full-scope simulation of the TN24P cask without using the 

concept of effective thermal conductivity or the porous media approach for the fuel region. The 

results from the FLUENT simulations were compared with experimental data on COBRA-SFS 

(Spent Fuel Storage) results. These comparisons showed good agreement with the observed 

deviation of the peak guide tube temperature, in that the temperature was within the uncertainty 

range of the thermocouples used in the experiment. However, there were deviations between 

the FLUENT and COBRA-SFS predicted velocity distributions. Cha et al. [2009] performed a 

preliminary thermal safety analysis of a commercial dry storage facility designed for the Korea 

standard SNF, in which they calculated the surface temperature of the storage canister 

corresponding to the SNF clad temperature.

Banken et al. [2010] used CFD analysis to validate the system designs for SNF storage and 

transportation applications. They explored the analytical steps necessary to achieve validation 

of the FLUENT CFD code models against the test results of Transnuclear’s SNF dry storage

systems and the CFD code model’s subsequent implementation in the analysis of another 

higher capacity system. The analysis found that the choices of meshing scheme, 
Introduction

1-3

temperature-dependent thermal properties, and turbulence modeling approach all influence 

the solution results. 

Adkins et al. [2013] presented thermal modeling studies for active storage modules in the 

Calvert Cliffs Nuclear Power Station Independent Spent Fuel Storage Installation (ISFSI). The 

thermal analysis was performed using the STAR-CCM+ commercial software, and the models 

developed for the specific ISFSI modules yielded temperature predictions in actual storage 

conditions for the concrete structure and the dry shield canister (DSC) and its contents, 

including preliminary estimates of fuel cladding temperatures for the SNF. The paper also 

presented a sensitivity analysis on boundary conditions to evaluate the effect of the modeling 

approach on predicted temperatures and temperature distributions. The results of this work 

demonstrated that existing CFD modeling tools can be used to obtain reasonable and accurate 

detailed representations of SNF storage systems with realistic decay heat loadings.

RES recently published NUREG–2152 [NRC 2012b], developed in cooperation with the Division 

of Spent Fuel Storage and Transportation (SFST) under NMSS. This NUREG provides 

practical advice and suggestions for reviewing CFD methods used in vendor applications of dry 

storage casks and for conducting high-fidelity CFD simulations of dry storage cask systems. It 

outlines procedures, analysis methods, and assumptions used for CFD simulations of dry 

storage cask systems. The procedures and the guidelines explained in NUREG–2152 for 

setting up the computational mesh and selecting the numerical schemes and turbulence models

were adopted in the current work.

1.3 Objective and Scope of Work

The objective of this work was to realistically model the thermal behavior of SNF in dry cask 

storage systems for periods of up to 300 years. The model was benchmarked with data from 

the Calvert Cliffs Nuclear Power Station ISFSI, and parametric studies and uncertainty analyses 

appropriate to the model were performed. The scope of this investigation consisted of three 

tasks. Task 1 was development of a three-dimensional (3-D) model of the NUHOMS® ventilated 

casks in the Calvert Cliffs Nuclear Power Station ISFSI. The model was benchmarked against 

measured profiles reported by a consortium of organizations including Pacific Northwest 

National Laboratory (PNNL), Oak Ridge National Laboratory (ORNL), Sandia National 

Laboratory (SNL), Electric Power Research Institute (EPRI), AREVA, and Constellation Energy 

[Suffield et al. 2012]. Task 2 consisted of simulating the thermal evolution of the casks for a 

period of 300 years. The thermal evolution was simulated for a range of heat generation rates 

from the spent fuel and ambient environmental conditions. Task 3 consisted of an analysis of 

the uncertainty in the simulation results using the GCI methodology described by ASME [2009].

The results of the analyses are presented and documented so that additional cases covering 

other dry cask storage system designs can be developed in a consistent manner for future 

licensing applications. 

1.4 Description of the Dry Cask Storage Modules

(HSM–1 and HSM–15) 

The dry cask storage modules at the Calvert Cliffs Nuclear Power Station ISFSI are a version of 

the NUHOMS design developed by Transnuclear, Inc., in which the spent fuel is sealed within a 

DSC that is loaded into a horizontal storage module (HSM). The DSC provides confinement, an 

inert environment with helium filling the void spaces, structural support, and criticality control for

PWR fuel assemblies [Banken et al. 2010]. The horizontal storage module is a reinforced 
Introduction

1-4

concrete storage overpack that provides for the passive removal of the spent fuel decay heat, 

environmental and seismic protection for the DSC, and radiological shielding from the DSC’s 

contents [Banken et al. 2010]. Figure 1-1 shows a typical NUHOMS HSM at the Calvert Cliffs 

Nuclear Power Station ISFSI. The figure shows the inlet vents and the front cover of the outlet 

vent structure. This figure also highlights the general arrangement of the storage modules on 

the ISFSI base pad. The specific storage modules considered in the present analysis 

are HSM–1 and HSM–15. Both HSM–1 and HSM–15 contain a DSC loaded with 24 CE 14 × 

14 assemblies. HSM–1 was loaded in November 1993, and the decay heat load for the DSC in 

HSM–1 calculated as of June 2012 was 4.1 kW [Suffield et al. 2012]. HSM–15 was loaded in 

November 1996, and the decay heat load for the DSC in HSM–15 was 10.8 kW at the time of 

loading and was calculated as 7.6 kW as of June 2012 [Suffield et al. 2012]. HSM–15 is one of 

the center units in a 2 × 6 array, while HSM–1 is a corner unit. 

Figure 1-2 shows the main components of the concrete structure, which has approximately a 

length of 6.1m [20.1ft], height of 6.5m [21.ft] and a width of 3.2m [10.6ft], (e.g., the air inlet, 

outlet, and plenum) and the position of the DSC shell. The ventilated storage system depends 

on natural convection for cooling. Cool air enters the system through the inlet vent and 

circulates in the plenum. Inside the concrete structure, the ventilation air absorbs decay heat, 

warming and expanding. The warmer, lighter air leaves the system through the outlet vent 

shown in Figure 1-2. To protect the concrete wall from the heat, a thin metallic heat shield is 

placed on the top and side internal walls of the concrete module. Internal components of 

the DSC shell, such as the spacer plates, tie rods, and fuel assemblies, are highlighted in 

Figure 1-3. The DSC is backfilled with pressurized helium to provide an inert atmosphere and 

enhance heat removal. Heat transfer from the fuel assemblies to the surrounding components 

and backfill helium takes place by conduction, convection, and radiation. The spacer plates 

shown in Figure 1-3 enhance the heat transfer and provide mechanical support to the fuel 

assemblies. The tie rods provide structural integrity and support to the system. The front and 

back covers of the DSC provide shielding. The DSC is placed on a pair of rails within the HSM 

module, as shown in Figures 1-2 and 1-3. 

Figure 1-1. NUHOMS horizontal storage modules in Calvert Cliffs Nuclear Power Station 

ISFSI [from Suffield et al. 2012]
Introduction

1-5

Figure 1-2. The concrete storage module and DSC shell with support rail for unit HSM–1

Figure 1-3. DSC shell and internal components: spacer plates, fuel assemblies, and tie 

rods for unit HSM–1

2-1

2. NUMERICAL MODELING METHODS

This chapter describes the mathematical and numerical models used to perform the simulations;

the rationale for selecting these models; assumptions, methods, and rationale used for mesh 

generation; and the basis for the material properties used in the simulations. Most of the 

numerical models used in the present work, the turbulence models used in the sensitivity 

analysis, the basis for mesh generation, and the material properties used in the simulations are 

based on guidelines provided in NUREG–2152, “Computational Fluid Dynamics Best Practice 

Guidelines for Dry Cask Applications” [NRC 2012b].

2.1 ANSYS-FLUENT Fluid Flow Solver 

The commercial computational fluid dynamics (CFD) package ANSYS-FLUENT® Version 14.5

[ANSYS 2012] was used for the simulations in the present study. ANSYS-FLUENT solves

the governing integral equations for the conservation of mass and momentum, and 

(when appropriate) for energy and other scalars such as turbulence and chemical species

concentrations. ANSYS-FLUENT uses a control-volume-based technique to convert a general 

scalar transport equation to an algebraic equation that is solved numerically. The following 

steps are used to solve the equations.

(i) Division of the domain into discrete control volumes using a computational grid.

(ii) Integration of the governing equations on the individual control volumes to construct 

algebraic equations for the discrete dependent variables (“unknowns”), such as 

velocities, pressure, temperature, and conserved scalars.

(iii) Linearization of the discretized equations and solution of the resultant linear equation 

system to yield updated values of the dependent variables.

2.2 Governing Equations 

For numerical simulation of all types of flows, ANSYS-FLUENT solves conservation equations 

for mass and momentum. For flows involving heat transfer or compressibility, an additional 

equation for energy conservation is solved. For flows involving species mixing or reactions, a 

species conservation equation is solved or, if the non-premixed combustion model is used, 

conservation equations for the mixture fraction and its variance are solved. Additional transport 

equations are also solved when the flow is turbulent. Only the governing equations relevant to 

gas flow and heat transport in the dry storage canister (DSC) are presented in this section.

2.2.1 Continuity and Momentum Equations

The equation for conservation of mass, or the continuity equation, can be written as follows: 

Equation (2-1)

Equation (2-1) is the general form of the mass conservation equation and is valid for 

incompressible as well as compressible flows. Here ρ is the density and is the velocity. The 

term defines the momentum. The source, Sm, is any user-defined source—for example, the 

mass added to the continuous phase from the dispersed second phase (due to vaporization of 

liquid droplets). 

( ) Sm t

+ ∇ • ρυ = ∂

∂ρ 

υ



ρυ


Numerical Modeling Methods

2-2

The conservation of momentum in an inertial (nonaccelerating) reference frame is described by

Equation (2-2)

where p is the static pressure, τ̿is the stress tensor, and ρg⃑ and F⃑ are the gravitational and 

external body forces, respectively. F⃑ also contains other model-dependent source terms, such 

as porous media and user-defined sources. 

The stress tensor τ̿is given by

Equation (2-3)

where ν is the molecular viscosity, I is the unit tensor, and the second term on the right-hand 

side is the effect of volume dilation. 

2.2.2 Energy Equation

The energy transport equation is given by

Equation (2-4)

where

= Unsteady term

= Convection term

= Conduction term

= Species diffusion term

= Viscous dissipation term

= Source term

= Velocity magnitude

In Equation (2-4), V is the velocity, keff is the effective thermal conductivity, h is the enthalpy, and 

J is the diffusion constant. 

The energy, E, per unit mass is defined as

Equation (2-5)

The numerics used in the present simulations, including the spatial scheme, the temporal 

scheme, and turbulence models, are explained in Section 2.5.

( ) ( ) p ( ) g F

t

     + ∇ • = −∇ + ∇ • + +

∂

∂

ρυ ρυυ τ ρ

( ) 



 



 = ∇ + ∇ − ∇ • I T τ ν υ υ υ   

3

2

( ) [ ( )] h

j

eff j j eff V E p k T h J V S

t

E

+



 





+ ∇ • + = ∇ • ∇ − ∑ + •

∂

∂

ρ τ ρ

( )

t

E

∂

∂ ρ

∇ •[V(ρE + p)]

keff ∇T

j j h J

τ eff •V

h S

V

2

2

p V E = h − +

ρ
Numerical Modeling Methods

2-3

2.3 Solid Modeling and Mesh Generation for

HSM–1 and HSM–15 

The two major differences between horizontal storage modules HSM–1 and HSM–15 in the 

Calvert Cliffs Independent Spent Fuel Storage Installation (ISFSI) are the heat load and the 

location of the module within the storage array. The HSM–15 module includes spent fuel with 

the highest decay heat load values in the entire Calvert Cliffs ISFSI [Adkins et al. 2013;

Suffield et al. 2012]. The DSC located within the HSM–1 module was the first canister loaded at 

the Calvert Cliffs ISFSI. According to Adkins et al. [2013], the decay heat load for the DSC in 

HSM–15 was 10.8 kW during the time of initial loading and the corresponding decay heat load 

for the DSC in HSM–1 was 5.8 kW (nearly 50 percent lower than HSM–15). HSM–15 is a 

central module in the 2 × 6 array of modules as shown in Figure 2-1, with only the top and front 

exposed to the atmosphere. The HSM–1 module is located on a corner of the 2 × 6 array and 

therefore has one side exposed to the atmosphere, in addition to the front and top. The 

exposed side wall of the HSM–1 module is much thicker than the interior walls between 

modules. External convection, insolation, and thermal radiation exchange with the environment 

were included in the model for the exposed side boundary of HSM–1.

The solid modeling and mesh generation for HSM–1 and HSM–15 include the mesh generated 

for the module internals, including the air inlet and air outlet vents; DSC support structures 

inside the concrete modules; and the heat shields (thin steel sheets) on the side walls and 

ceilings to protect the concrete walls from the thermal load of the DSC. The generated model 

and computational grid also include the DSC interior structure, comprising the spacer disks, tie

Figure 2-1. Aerial view of the Calvert Cliffs Nuclear Power Station ISFSI housing HSM–1 

and HSM–15 modules [imagery© U.S. Geological Survey, Map Data© Google 

Inc. 2013]
Numerical Modeling Methods

2-4

rods, fuel assemblies, and the DSC top and bottom end caps. The stainless steel guide sleeves 

surrounding the fuel, along with the air gap between the sleeve and fuel assemblies, were also 

explicitly modeled. Figure 2-2 illustrates the geometry of the DSC and the three-dimensional 

(3-D) model of the DSC and its components. Figure 2-3 shows a 3-D schematic view of 

HSM–1. The construction of HSM–15 (not shown) is very similar to HSM–1, except that both 

side walls are the same thickness. Figure 2-4 shows a 3-D view of the solid model developed 

for DSC internals, such as the support rod, guide sleeve, and spacer disk, and the supporting 

structure, such as guide rails. Figure 2-5 shows the solid model developed for the fuel 

assemblies, tie rods, and the spacer plates.

2.3.1 Grid Generation

The same basic grid generation methodology was used for both HSM-1 and HSM-15. The 

computational domain consisted of the DSC geometry, the HSM module, the airflow path, the 

inlet and outlet vent, and a portion of the surrounding ambient environment. The baseline mesh 

consisted of 3.5 million cells over the entire computational domain. The computational grid was 

generated following guidance in NUREG–2152 [NRC 2012b]. To develop a controlled and 

systematic grid, the mesh density along the edges of the geometrical entities and the mesh 

refinement parameters such as stretching ratios were specified. Hexahedral grids were utilized 

to mesh the fluid and solid volumes in the domain. In most cases, the volume meshing was 

done by sweeping a two-dimensional (2-D) quadrilateral mesh along a specific direction. 

Figure 2-6 shows such a 2-D meshed surface that was used to generate the grid for the internal 

components within the DSC shell. The mesh generation method used in this study created a 

Figure 2-2. Schematics showing the geometry of the DSC [courtesy AREVA, PNNL]
Numerical Modeling Methods

2-5

Figure 2-3. Three-dimensional schematic view of the HSM–1 module

Figure 2-4. Three-dimensional view of the DSC internal showing the support rod, guide 

sleeve, spacer disk, and guide rails
Numerical Modeling Methods

2-6

Figure 2-5. Three-dimensional schematic for the fuel assemblies, tie rods, and spacer 

plates within the DSC shell

Figure 2-6. Two-dimensional cross-sectional view showing the computational grid for 

the fuel region
Numerical Modeling Methods

2-7

compact mesh that was amenable to systematic and controlled refinement, which was needed 

for the grid convergence index study described in Chapter 5. 

The mesh spacing within the DSC region and in the region of the airflow path around the DSC 

was resolved more finely for accurately (i) simulating the airflow through the airflow path around 

the DSC and (ii) capturing the temperature and velocity gradients in the DSC region. A 

fine-resolution mesh was used in the near-wall region to facilitate the use of the low Reynolds 

number k-ε model. Ideally, a y+ of 1 should be used to resolve velocity and temperature 

gradients up to the wall without the use of any wall functions. From the baseline simulations, it 

was found that the surface average wall y+ for the DSC shell outer wall was approximately 2.97. 

Subsequently, a wall mesh refinement study was performed to understand the effect of near 

wall mesh size on the computational results (Section 3.3.1). It was found that near wall mesh 

refinement did not significantly affect the results and the baseline mesh was deemed 

appropriate for use in the study. In the regions where the velocity and temperature gradients 

are mild (such as in the concrete walls and in the airflow regions far from the DSC), a relatively 

coarse mesh was used for computational efficiency.

Figures 2-7 through 2-10 show the exterior and interior meshes for the HSM–1 module. The 

schematic for HSM–15 is the same. Figure 2-7 shows the 3-D exterior view of the 

computational grid for the HSM–1 module. Figure 2-8 shows a transverse cross-sectional plane

of the domain and highlights the internal mesh used for the air volume within the concrete 

structure surrounding the DSC shell. Figure 2-9 shows a cross-sectional plane highlighting the 

general pattern of the internal mesh used in the axial plane. Figure 2-10 provides a closer look 

at the axial mesh generated within and around the DSC shell, fuel assemblies, spacer plates, 

and helium backfill region. The interior mesh, including the DSC within the HSM–1 module, is 

illustrated in Figures 2-11 and 2-12. Figure 2-11 shows the mesh generated on the DSC shell 

and support structure, while Figure 2-12 illustrates the mesh generated for the fuel assemblies, 

spacer plates, and tie rods. 

The computational grid for HSM–15, with its thinner side wall, consisted of 3.2 million cells. The 

grid generation for the HSM–15 module followed the same principles that were used to generate 

the computational grid for the HSM–1 module and is not discussed in detail. As a part of the 

numerical uncertainty analysis, a detailed grid convergence index (GCI) study was performed 

(Chapter 5). The GCI study requires at least four mesh levels that are developed based on 

systematic refinement or coarsening of a baseline mesh. The existing mesh for HSM–1 and 

HSM–15 was used as the baseline mesh. Two levels of refinement and one level of coarsening 

were performed to create the four meshes necessary for the GCI study. The approach for mesh 

refinement and coarsening used in this study was adopted from the guidelines for best practice 

use of CFD for dry storage cask analysis provided in NUREG–2152 [NRC 2012b]. The total 

number of cells for each mesh level is shown in Table 2-1.

2.4 Boundary and Operating Conditions

The following subsections describe the boundary conditions used for HSM–1 and HSM–15 in 

the current simulations. Though mostly similar, some of the boundary conditions for HSM–1 

and HSM–15 were applied differently because of the different locations of the modules. As 

mentioned in an earlier section, the HSM–15 module is located in the middle of the 2 × 6 array 

of modules and, hence, an adiabatic wall (no heat flow) boundary condition was applied for both 

side walls assuming similar decay heat of the casks at both sides. The HSM–1 module is 

located at one end (edge) of the 2 × 6 array of modules and therefore has one side exposed to 

the ambient condition and the other is assumed to be symmetric. In reality, some heat transfers 
Numerical Modeling Methods

2-8

Figure 2-7. Three-dimensional view of the computational grid (exterior view) of the 

HSM–1 module

Figure 2-8. Transverse slice through center line showing computational grid of HSM–1, 

highlighting the airflow path
Numerical Modeling Methods

2-9

Figure 2-9. Planar slice through mid-section of the DSC showing computational mesh of 

HSM–1 module in transverse direction

Figure 2-10. Enlarged view of the computational grid in the guide sleeve, support rod,

and spacer plate region
Numerical Modeling Methods

2-10

Figure 2-11. Three-dimensional view of the computational grid for the DSC internal and 

guide rails

Figure 2-12. Three-dimensional view of the DSC internal showing the fuel region, 

support rod, guide sleeve, spacer disk, and guide rails
Numerical Modeling Methods

2-11

Table 2-1. Total number of cells for different grid levels used in the GCI study

Total number of cells (in millions)

Refinement/Coarsening HSM–1 HSM–15

Refinement Level 2 24.51 22.63

Refinement Level 1 9.21 8.51

Baseline 3.52 3.27

Coarsening Level 1 1.22 1.13

from the central module to the adjacent corner modules through the side walls, depending upon 

the heat load of adjacent modules.

Once the model development and mesh generation for all the components were complete, fluid 

and solid regions were specified within the storage module and the canister. An air region was 

specified within the concrete structure of the storage module that surrounds the DSC. A 

helium-filled region was specified within the DSC. The solid region was composed of the 

concrete structure, DSC shell, support structure (e.g., rails), and DSC shell internals 

(e.g., spacer plates, fuel sleeves, tie rods, and end caps). The material for each of the solid 

regions was also specified. Suitable boundary conditions were specified at the exposed outer 

surface of the domain. An ambient environment was assumed to exist outside the HSM.

On the exposed boundary of the concrete (for both HSM–1 and HSM–15), a combination of 

external convection and thermal radiation exchange with the environment was allowed. The 

heat transfer factors for the convection heat transfer were calculated based on the correlation 

for free convection from the vertical flat plates and horizontal surfaces to the surrounding still air. 

The natural convection correlations used to calculate the heat transfer factors are presented in 

Appendix A. Radiation from the side surfaces and top walls was also included and an emissivity 

of 0.9 for the concrete external surface was specified. An ambient temperature was specified at 

301 K [82 °F] for all the simulations, representing equilibration to late June air temperatures at 

the site to represent conditions during the temperature measurements. On the bottom surface, 

which is the interface between the module and underlying soil, conduction through the base to a 

concrete pad and underlying soil was specified by defining appropriate thermal resistance. 

This thermal resistance is determined as an equivalent conduction through 1.8 m [6 ft] of 

soil to a temperature of 288 K [58.73 °F], with a typical soil conductivity of 0.52 W/m-K 

[8.34 BTU/ft-sec-°F] [Zigh and Solis 2008]. The temperature at depth within the soil is the mean 

annual air temperature at the site. Solar heat loads were applied at the outer exposed surfaces 

of the modules. The regulations for packaging and transportation of radioactive material 

(10 CFR Part 71, Section 71.71) specify, for a 12-hour period, a solar load of 800 g-cal/cm2

[2,949.35 BTU/ft2

] for horizontal surfaces, 200 g-cal/cm2 [737.337 BTU/ft2] for flat non-horizontal 

surfaces, and 400 g-cal/cm2 [1,474.675 BTU/ft2

] for curved surfaces. Because there are no 

curved external surfaces, the solar load specification used the values for horizontal and vertical 

flat surfaces. A heat generation rate that produces a heat flux equivalent to the required solar 

loading was specified on the external surfaces in the ANSYS-FLUENT model.

The storage modules at Calvert Cliffs have a single air inlet on the bottom of the front face and 

two exit vents at the top of the module. These vents are protected by a number of screens, 

which provide resistance to airflow and therefore cause a pressure drop. To appropriately 

account for the pressure drop, the inlet and outlets are modeled as vents in ANSYS-FLUENT, 

where the pressure loss factor is specified as a piecewise linear function of normal velocity. 

The terms of this equation were derived based on the pressure drop equation for 

incompressible flows through screens [Perry and Green 1997]. This method is described in 

Appendix A. 
Numerical Modeling Methods

2-12

The back surface of the storage module for both HSM–1 and HSM–15 was assumed to be

symmetric to account for the presence of another storage module with similar decay heat. Both 

the side surfaces of HSM–15 and one of the side surfaces of HSM–1 were also specified to 

have a symmetric boundary. Though specification of symmetry implicitly assumes that the 

thermal loads of the adjoining modules are the same, this assumption enables estimating the 

temperature field with reasonable computing resources and without including all storage 

modules in the computational domain.

The ambient pressure was specified as 101,325 N/m2 [2,116 lbf/ft2

]. The operating ambient air 

density was specified as 1.173 kg/m3 [0.072 lbm/ft3

], which was calculated using the ideal gas 

law based on the ambient pressure and temperature. For the internal surfaces that are not 

exposed to the ambient environment, coupled boundary conditions were used at the solid–fluid 

interface. Based on the material of the solid body, emissivity values were specified for the 

surface. The emissivity values associated with specific materials used in the model are listed in 

Table B–7 of Appendix B. These values were obtained from engineering handbooks and 

textbooks, and were not directly measured or tested for the actual cask configuration.

2.5 Numerics, Radiation Heat Transfer, and Turbulence Models for 

the Simulations

The pressure-based, steady-state solver of ANSYS-FLUENT was used in simulating the flow 

field of the HSM–1 and HSM–15 modules. The velocity coupling was done using the SIMPLE 

algorithm. For spatial discretization, second-order upwind spatial differencing was used for all 

variables except the pressure equation (continuity equation), where a body force weighting 

factor was used. Gradients were calculated with a Green-Gauss node-based method. Suitable 

under-relaxation factors were used for all the variables to optimize numerical stability and 

achieve faster convergence.

Thermal radiation between surfaces within the HSM module elements and DSC components is 

vital to the overall heat transfer process. In the present study, the discrete ordinate model (DO) 

was used for thermal radiation modeling. In this approach, the radiative transfer equation (RTE) 

for an absorbing, emitting, and scattering medium is solved for a finite number of discrete solid

angles. The discrete ordinate model transforms the RTE into a transport equation for the 

radiation intensity in the spatial coordinates. The number of equations the DO model solves is 

controlled by the angular discretization. To solve the HSM–1 and HSM–15 problems, four 

angular discretizations were used in each direction of the spherical coordinate system [theta (θ) 

and phi (φ)]. The number of energy iterations per radiation iteration was limited to 10 for optimal 

convergence. The modeled storage systems do not have any semi-transparent media; hence 

explicit DO-energy coupling was not needed to enhance convergence.

The Reynolds number for the airflow within the concrete module was estimated using the 

average air velocity and density from initial HSM–1 runs and a characteristic length representing 

the flow path. The calculated Reynolds number was approximately 8,399, which is significantly 

higher than the threshold of 4,000 for producing turbulent internal flows. As the physics of the 

problem involves natural convection, the Rayleigh number of the flow was determined based on

the DSC shell diameter, average temperature of the air and DSC shell, and average density and 

fluid properties. The calculated Rayleigh number was 8.45 × 109

, which is higher than the 

critical value for natural convection on vertical plates (109

). Because the flow within the storage 

module is different from simple vertical flat plate flows, the added complexity is expected to 

trigger an early transition to turbulence, thus the flow will be in the transitional turbulent regime 

for HSM–1. HSM–15 contains hotter fuels, so the flow velocity for HSM–15 will be higher and 
Numerical Modeling Methods

2-13

result in a larger Reynolds number. The higher HSM–15 temperature will yield a larger 

Rayleigh number, so the airflow within HSM–15 will also be in the transitional or turbulent 

regime. An appropriate turbulence model is needed to accurately predict local velocities and 

temperatures in the air stream. The low Reynolds number k-ε model with full buoyancy effects 

was used for the majority of the simulations performed in this study. In addition, a sensitivity 

study was performed to assess the effect of turbulence models on calculated quantities

(Section 3.3.2). The standard k-ε model with scalable wall function and buoyancy effects, SST 

k-ω model with low Reynolds number corrections, and SST k-ω model with buoyancy effects 

were used in this sensitivity study. The Reynolds number for the helium flow within the DSC 

shell was in the laminar region; thus, the helium flow was treated as a laminar zone in 

the simulations.

2.6 Fuel Assembly Model

HSM–1 and HSM–15 both contain 24 pressurized water reactor (PWR) fuel assemblies, each 

carrying an arrangement of CE 14 × 14 fuel rods. Each fuel assembly inside the DSC shell was

modeled as a porous medium using guidelines provided in NUREG–2152 [NRC 2012b] and as 

described in Appendix C. The fuel assembly model development consisted of three steps: (i) a 

porous media model calculation to determine flow resistances, (ii) an effective thermal 

conductivity calculation, and (iii) a fuel assembly decay heat loading specification. Each of 

these steps is detailed in separate appendices. 

2.6.1 Porous Media Calculation

The porous media model was used in the HSM–1 and HSM–15 models to represent flow inside 

the fuel assemblies, which are each treated as a continuous porous medium with predetermined 

resistance specified in the axial and radial flow directions. These parameters are used as 

surrogates for the pressure loss caused by the presence of fuel rods. In essence, the porous 

media model adds a momentum sink in the governing momentum equations. This 

approximation significantly reduces the computational mesh requirements needed to explicitly 

model each fuel rod in the assembly. 

A frictional resistance and an inertial resistance factor must be specified for the porous media 

model in ANSYS-FLUENT. 3-D CFD simulations of an isolated fuel assembly were performed

to determine these resistance factors. The 3-D simulations model the complete length of a fuel 

assembly containing CE 14 × 14 fuel rods, represented as solid cylindrical objects. Simulations 

were performed for a number of inlet velocities encompassing the range of velocities expected. 

From each simulation, the pressure drop across the fuel assembly was determined, then 

frictional and inertial resistances were calculated from the set of calculated pressure drops and 

input velocities. The porous media model calculations are described in Appendix C.

2.6.2 Effective Thermal Conductivity Calculation

The fuel assemblies were modeled as porous regions, assumed to have anisotropic thermal 

conductivity. Based on the fuel arrangement within an assembly, the effective radial and axial 

thermal conductivities were determined in a separate, stand-alone calculation. 

The temperature-dependent effective thermal conductivity was determined in the axial direction 

using a simple weighted-area-average method. The fuel (UO2) material, helium fill gas, and 

Zircaloy cladding of the fuel rod were considered in the averaging process. In the radial 

direction, the effective thermal conductivity method [Bahney and Lotz 1996] was used to obtain 

the conductivity of the fuel as a function of temperature. To calculate the 
Numerical Modeling Methods

2-14

temperature-dependent thermal conductivity, a separate 2-D finite volume model, using 

ANSYS-FLUENT 14.5, was built to represent a planar section of the fuel assembly. This 2-D

model considered each fuel rod along with fuel cladding and the backfill helium gas. Every fuel 

rod in the assembly was assumed to have the same heat load. A heat load value 

representative of the conditions expected in the present study was used in the analysis. A 

series of simulations was performed to determine the temperature-dependent effective radial 

thermal conductivity of the fuel assembly. The calculated axial and radial effective thermal 

conductivity values were used in ANSYS-FLUENT as temperature-dependent conductivity 

functions using an orthotropic distribution. The effective thermal conductivity calculations are 

described in Appendix D.

2.6.3 Fuel Assembly Loading

The decay heat load was defined for each of the 24 active fuel regions of the CE 14 × 14 

assemblies within the DSC for both HSM–1 and HSM–15 based on the decay heat load data 

from the Calvert Cliffs ISFSI provided by NRC. The fuel structure and fuel identification location 

for HSM–1 and HSM–15 are shown in Figures 2-13 and 2-14, respectively. The decay heat 

load was applied as a uniform volumetric heat generation rate over the active fuel length within 

each fuel region. This volumetric heat generation rate was modified to include the axial 

variation by multiplying the volumetric heat generation rate by an axial peaking factor. The axial 

variation of heat load was implemented as a user-defined function in the ANSYS-FLUENT 

solver. The decay heat calculation and fuel assembly loading for the HSM–1 and HSM–15 

modules are detailed in Appendix E. 

2.7 Material Properties Used in the Simulations

In general, temperature-dependent material properties were used in the simulations. Concrete, 

stainless steel, carbon steel, and lead are the solid materials used in the model, and air and 

helium are the fluid materials. The material properties are detailed in Tables B–1 through B–7 

of Appendix B.
Numerical Modeling Methods

2-15

Figure 2-13. Fuel assembly loading identification structure for HSM–1 module

Figure 2-14. Fuel assembly loading identification Structure for HSM–15 module

3-1

3. BENCHMARK MODELING AND SENSITIVITY ANALYSES

This chapter discusses the results of the benchmark simulations for horizontal storage modules

HSM–1 and HSM–15. The computed temperature data are compared with field measurements. 

This chapter also presents results from a series of sensitivity studies on the HSM–1 simulation, 

performed to assess the effects of turbulence models, porous media resistance, and insolation. 

An additional simulation was performed using an extended computational domain to include the 

surrounding atmosphere to verify the convective heat transfer coefficient used at the outer walls 

in the baseline study.

3.1 Description of Field Temperature Measurement 

The field temperature measurements for the HSM–1 and HSM–15 modules in the Calvert Cliffs 

Nuclear Power Station Independent Spent Fuel Storage Installation (ISFSI) were taken on 

June 27 and 28, 2012. The direct temperature measurements were conducted under the Used 

Fuel Disposition Campaign of the U.S. Department of Energy (DOE) [Suffield et al. 2012]. 

Visual inspections and surface sampling were also performed. Due to safety constraints and 

physical constraints related to accessing the different regions of the DSC, reliable temperature 

measurements were obtained only on the exposed face of the canister base and a short 

distance along the canister side [Suffield et al. 2012]. The front door of each module was 

opened to obtain the temperature measurements. The HSM–15 module was open for about 

20 minutes, with a single temperature measurement obtained on the exposed base of the DSC

[Suffield et al. 2012]. HSM–1 was open for approximately 160 minutes, with 16 temperature 

measurements obtained on the base and sides of the canister [Suffield et al. 2012]. 

The temperatures were measured with an Omega “all-in-one” hand-held thermometer, 

Model 450-AET with Type E Chromega®-constantan thermocouple, or a Model 450-APT 

platinum resistance temperature detector (RTD) thermocouple [Suffield et al. 2012]. The 

observation points and comparisons between the measured and simulated data are presented 

in the following sections.

3.2 Results of Baseline Simulations

Computed results from the baseline steady-state simulations of HSM–1 and HSM–15 are 

presented in this section. The decay heat loads for each fuel assembly used in the baseline 

steady-state simulations are listed in Appendix E. The ambient temperature considered in the 

simulations is 27.9 °C [301 K, 82.19 °F], corresponding to the mean air temperature for the day 

when measurements were taken. 

The spatial distributions of simulation results are presented using three-dimensional (3-D) plots 

and cross sections. Figure 3-1 displays the axial and transverse cross sections used to display 

temperature and velocity patterns for both the HSM–1 and HSM–15 simulations. Both of these 

cross sections pass through the centroid of the DSC.

3.2.1 HSM–1 

The total decay heat load for the dry storage canister (DSC) in the HSM–1 module was 

calculated at approximately 4.1 kW [13,990 BTU/hr]. This decay heat was used along with the 

volume of each fuel assembly region to determine the volumetric heat generation rate input in 

the simulations. This volumetric heat generation rate was applied over the active fuel length 

within each fuel region and was modified to include the axial variation through an axial peaking 
Benchmark Modeling and Sensitivity Analyses

3-2

Figure 3-1. Axial and transverse cross-sectional planes used to illustrate temperature 

and velocity distributions

factor. The peaking factor is constant at a maximum value over the central 60 percent of the 

axial length, and drops approximately linearly to the end of the canister to roughly half of the 

peak value. The peaking factor calculations are described in Appendix E.

3.2.1.1 Temperature and Velocity Patterns

Figures 3-2 and 3-3 present simulated temperatures on the DSC shell surface and on internal 

metal surfaces within the DSC for the HSM–1 module. Figure 3-2 shows that temperatures 

peak at approximately 75 K [135 °F] above the ambient temperature in the upper middle section 

of the DSC, and drop to approximately 15 K [27 °F] and 20 K [36 °F] above ambient at the back 

and front ends, respectively. The back end of the DSC is cooler than the front end because of 

enhanced ventilation flow and the location of the active fuel zone. The high-temperature region 

is located at the upper middle section of the DSC because (i) decay heat loads are largest in the 

central portions of the individual fuel assemblies and (ii) temperature-driven natural convection 

flow of the helium coolant inside the DSC transports heat up from the core of the DSC. The 

helium coolant forms two natural convection cells, with warming helium rising in the center of 

the DSC and cooling helium falling along the DSC side walls as heat is transferred through the 

DSC wall to the ventilation air. It can also be observed that the support rails, at only 10 K 

[18 °F] above ambient, are relatively cool compared to the other parts of the DSC shell, losing 

heat by conduction through contact with the concrete side walls.

Figure 3-3 illustrates the temperature distributions for the fuel region, spacer plates, and support 

rods. Temperature is highest at the upper part of the fuel region and peaks at 85K [153 °F] 

above ambient, which is approximately 10K [18 °F] higher than the peak DSC shell temperature. 
Benchmark Modeling and Sensitivity Analyses

3-3

Figure 3-2. Contours of DSC shell surface temperature for the HSM–1 baseline case

Figure 3-3. Temperature distributions for the fuel region, spacer plates, and support rods

for HSM–1 baseline case
Benchmark Modeling and Sensitivity Analyses

3-4

This can again be attributed to buoyancy driven natural convection transporting warmer helium 

upwards, resulting in a high temperature region at the upper section of the fuel region. The high 

temperature region is also located near the middle of the fuel region in the axial direction, with 

cooler temperatures at the front and back of the assemblies. As observed in connection with 

Figure 3-2, this pattern results from the peaking factor and heat load distribution along the 

axial direction.

Figures 3-4 through 3-7 illustrate the spatial distribution of temperature in (i) fluids in the 

transverse cross section, (ii) solids in the transverse cross section, (iii) fluids in the axial cross 

section, and (iv) solids in the axial cross section, respectively. The warm temperatures in the 

upper middle DSC are clearly evident in all four figures. Figure 3-4 shows that the ventilation air 

is approximately 10 K [18 °F] warmer in the poorly connected volume between the HSM–1 door 

and the front end of the DSC than in the well-connected volume at the back end of the DSC, 

and shows that the ventilation air increases in temperature by 45 K [81 °F] from the inlet to 

the outlet. A comparison of Figures 3-6 and 3-7 illustrates the sharp drop in temperature 

{30 K [54 °F] to 50 K [90 °F]} across the DSC shell.

Figures 3-8 through 3-10 show the predicted velocity field of the HSM–1 module, using airflow 

streamlines in Figure 3-8 and air and helium velocity vectors on the transverse and axial cross 

sections in Figures 3-9 and 3-10.

The airflow streamlines in Figure 3-8 illustrate the ventilation process for the dry storage system. 

Cool air enters the system at the inlet and gains heat from the DSC. The air expands as it 

warms, with the less-dense air inducing flow out the top vents and causing air to flow within the 

cavities around the DSC. The relatively large velocities in the cavity along the DSC walls and 

above the top of the DSC are examples of buoyancy-induced convection due to the strong 

heating from the DSC. The colored streamlines and velocity magnitude plots indicate that the 

fastest air movement is near the inlet and outlet, and in the narrow cavity below the DSC. Low 

velocities occur in the open region below the DSC, where the inlet jet creates a large 

recirculation zone. The poorly connected volume at the front end of the DSC has very low 

velocities and no streamlines enter the volume; the well-connected volume at the back end of 

the DSC also has relatively low velocities, but some streamlines do enter this zone. In general, 

the air flow within the concrete module and around the DSC shell is asymmetric as a result of 

design features such as the locations of the air inlet vents, outlet paths, and internal 

obstructions. This flow asymmetry contributes to the uneven cooling of the DSC shell and 

internal components in the axial direction, as illustrated previously in Figures 3-2 and 3-3.

Figures 3-9 and 3-10 show that the movement of helium within the DSC is highly localized into 

channels between the fuel assemblies and is roughly proportional to the heat loading. Within 

these channels, helium velocities are comparable to the ventilation air velocities outside the 

DSC shell, but helium is essentially stagnant within each fuel assembly. Warming helium 

moves upward within the channels and cooling helium drops along the DSC side walls.

3.2.1.2 Temperature Line Probes

The simulated temperature distribution on the DSC shell, heat shield, and fuel region can be 

examined by extracting values along predefined axial probe lines from the ANSYS-FLUENT 

results. Probe lines were selected to demonstrate relatively high and low temperatures on the 

DSC shell, on the heat shield, and within particular fuel assemblies. Figure 3-11 shows top and 

right probe lines for the DSC shell and heat shield, representing high and low temperatures, 

respectively. Figure 3-12 shows fuel assembly probe lines for relatively hot (fuel assembly 1) 
Benchmark Modeling and Sensitivity Analyses

3-5

Figure 3-4. Air and helium temperature distributions along the transverse cross section 

for the HSM–1 baseline case

Figure 3-5. Temperature distribution for the solid region at transverse cross section for 

HSM–1 baseline case
Benchmark Modeling and Sensitivity Analyses

3-6

Figure 3-6. Temperature distribution for the fluid region along the axial cross-sectional 

plane for the HSM–1 baseline case

Figure 3-7. Temperature distribution for the solid region along the axial cross-sectional 

plane for the HSM–1 baseline case
Benchmark Modeling and Sensitivity Analyses

3-7

Figure 3-8. Streamlines showing the airflow path and airflow circulation within the 

HSM–1 module

Figure 3-9. Velocity vectors at the mid-plane in the transverse cross section for HSM–1 

module
Benchmark Modeling and Sensitivity Analyses

3-8

Figure 3-10. Velocity vectors at the axial cross section for HSM–1 module

Figure 3-11. Line probe locations for extracting temperature distribution along the axial 

direction: locations shown on the DSC shell and on the heat shield
Benchmark Modeling and Sensitivity Analyses

3-9

Figure 3-12. Line probe locations for extracting temperature distribution along the axial 

direction in the fuel region: locations shown for fuel assemblies 1 and 23

and cool (fuel assembly 23) assemblies. The corresponding temperature profiles are shown in 

Figures 3-13 through 3-15.

Comparing Figures 3-13 and 3-14, it can be seen that the temperature difference between the 

top and right probe lines is essentially the same for the DSC shell and the heat shield. The right 

side of the DSC shell shows a fluctuating pattern with spikes (Figure 3-13), which can be 

attributed to the presence of spacer plates. In general, the spacer plates have a higher 

temperature than the surrounding helium. Hence, the proximity of spacer plates to the DSC 

shell wall enhances radiative heat transfer, resulting in localized high temperature zones. The 

relatively hot helium at the top of the DSC shell provides a more uniform temperature field next

to the shell, damping the peaks. In both the DSC shell and the heat shield, the temperature 

distribution is skewed towards the module door as a result of asymmetric ventilation air patterns 

and the location of the active fuel zone (see Figure 3-8).

Figure 3-15 contrasts temperatures along fuel assemblies 1 and 23. The assemblies are both 

hotter and more symmetric than the shell and heat shield. Assembly 1 is much hotter than 

assembly 23 because of natural convection within the DSC. The peak temperatures in these 

assemblies are not located at the center, but are biased towards the front of the module. This 

effect, more prominent in assembly 23, is primarily due to the distribution of heat load along 

the axial direction of the fuel assemblies (illustrated by the peaking factor distribution in 

Appendix E). However, this may also be related to uneven cooling due to asymmetric flow 

around the DSC shell, as highlighted in Figure 3-8.
Benchmark Modeling and Sensitivity Analyses

3-10

Figure 3-13. Computed temperature distribution along the axial distance for top of DSC 

shell and right of DSC shell

Figure 3-14. Computed temperature distributions along the axial distance for top of heat 

shield and right of heat shield
Benchmark Modeling and Sensitivity Analyses

3-11

Figure 3-15. Computed temperature distributions along the axial distance for fuel region: 

distribution shown for fuel assemblies 1 and 23 for HSM–1 baseline case

3.2.1.3 Comparison With Measured Data 

Table 3-1 presents a comparison of the computed temperature with measured temperatures at 

specified locations. The measurement locations were constrained by the length of the 

thermocouple probe, which reached a maximum of 1.016 m [3.33 ft] along the sides of the DSC.

Measurements were obtained at the top, sides, and support rail locations (Figures 3-16 and 

3-17). The temperature measurements were taken at three axial locations:

• 1.02 m [3.34 ft] from the exposed surface of the DSC surface,

• 0.51 m [1.67 ft] with the probe inserted half the length, and

• 0 m [0.0 ft].

From the description of the measurement procedure [Suffield et al 2012], it is apparent that the 

measured temperatures have significant uncertainty. Before the measurements, the front door 

of HSM-1 was open for approximately 160 minutes. This allowed hot air from within the module 

to escape and lower the temperature within the system. Therefore, the measured temperatures

would be lower than the temperatures of a normally operated closed module. Temperature

measurements were taken by touching a handheld device to the canister surface [Suffield et al 

2012]. It is not clear if the probe was in contact with the surface for a period long enough to 

establish thermal equilibrium between the proble and surface. Without thermal equilibrium, the 

probe would record a temperature lower than the actual value. A review of measured data from 

Table 3-1 shows a relatively flat temperature distribution along the axial direction, especially at

the DSC shell top wall. Given the magnitude of the heat load and the axial distribution of the 

peaking factor, a flat axial temperature distribution is unlikely. This may indicate that, due to the
Benchmark Modeling and Sensitivity Analyses

3-12

Table 3-1. Comparison of measured and predicted temperatures for the baseline case 

simulations for HSM–1 module

Location 1 m [3.28 ft]

Measured value (K) 

[°F = 1.8 K-460]

Baseline case simulation 

value (K) [with insolation, 

low Re standard k-ε model 

with buoyancy] (K)

[°F = 1.8 K-460]

Under Grapple ring 317.44 324.18

0-m-side-right 315.22 323.04

0-m-side-top 319.11 325.18

0-m-side-left 313 322.68

0-m-rail-right 313.55 319.05

0-m-rail-left 314.11 319.05

0.51-m-side-right 315.77 325.74

0.51-m-side-top 320.22 336.94

0.51-m-side-left 313.55 324.94

0.51-m-rail-right 314.11 318.81

0.51-m-rail-left 314.66 318.7

1.02-m-side-right 315.22 334.35

1.02-m-side-top 321.33 346.89

1.02-m-side-left 315.22 333.39

1.02-m-rail-right 314.11 322.04

1.02-m-rail-left 315.22 321.97

Figure 3-16. Thermocouple locations shown at 0.0-, 0.51- and 1.02-m axial distance
Benchmark Modeling and Sensitivity Analyses

3-13

Figure 3-17. Thermocouple locations shown at the top, right side, left side, rail right, and 

rail left

difficulty of inserting the proble in the narrow space between the DSC shell and module 

concrete wall, complete contact between the probe and DSC shell wall was not established and 

the probe reported the air temperature rather than the solid surface temperature. Finally, the 

temperature measurement was performed only once. Without repeated measurements, the 

uncertainty in the measurement cannot be evaluated.

At each axial location, temperature measurements were taken at five different locations: at two 

horizontal sides and the top of the DSC shell, and at the two sides of the support rails. These 

locations are highlighted in Figure 3-17. In addition, a measurement was taken under the 

grapple ring, which is located on the cover of the DSC (Figure 3-16).

Table 3-1 illustrates that the simulated temperatures, in general, overpredict the measured 

temperatures. The minimum deviation between the computed and measured data is 

approximately 4 K [7.3 °F] and occurs at the middle section on the rail. The maximum deviation 

is around 25.5 K [46.0 °F] and occurs at the top interior section of the DSC shell. The difference 

between measured and observed temperatures increases with distance from the opening. As 

discussed previously in this section, there is significant uncertainty associated with the

measured temperatures and whether they are appropriate for comparison with the current 

predicted results. The simulated temperatures are in reasonable agreement with the measured 

temperatures near the cap of the DSC {i.e., at an axial location of 0.0 m [0.0 ft]}, but the 

temperatures measured at 0.51 and 1.02 m [1.67 and 3.34 ft] at the surface are significantly 

lower than the simulated temperatures. In addition, as discussed previously in this section, the 

thermocouple probe may not have made good contact with the DSC surface on the side of the 

DSC, in which case the measured temperatures at 0.51 m [1.67 ft] and 1.02 m [3.34 ft] may 

actually be air temperatures and not DSC surface temperatures. Suffield et al. [2012] 

expressed similar observations and views.
Benchmark Modeling and Sensitivity Analyses

3-14

3.2.2 HSM–15 

Computed results from the baseline simulation of the HSM–15 module are presented in this 

section. The baseline simulations were performed using an ambient temperature of 28 °C 

[301 K, 82.19 °F] and decay heat load values calculated for June 2012. The total heat load for 

HSM–15 in June 2012 was 7.6 kW [25,933 BTU/hr], compared to 4.1 kW [13,990 BTU/hr] for 

HSM–1. The decay heat load distribution for each fuel assembly region in the HSM–15 module 

is provided in Appendix E, using the same axial peaking factor distribution as was used for 

HSM–1. This decay heat was divided by the volume of each fuel assembly region and used as 

the volumetric heat generation rate input in the simulations. The HSM–15 module is located at 

the middle of the 2 × 6 array of modules, thus only the top and front are exposed to the 

atmosphere. In contrast, the HSM–1 module is located at one corner of the 2 × 6 array of 

modules and therefore has a side exposed to the atmosphere. Because the general patterns of 

temperature and flow are similar to those reported for HSM–1, only the differences between the 

HSM–15 and HSM–1 simulations will be discussed.

3.2.2.1 Temperature Field

Figures 3-18 through 3-26 present the simulated temperature field of the HSM–15 module. 

Figures 3-18 and 3-19 show the temperature distribution along the DSC shell surface and on 

internal surfaces, respectively. The highest simulated temperature for the HSM–15 DSC shell 

{391.46 K [245.2 °F]} is significantly higher than the highest temperature simulated on the DSC 

shell for the HSM–1 module {355.27 K [179.82 °F]}. This observation is consistent with the 

larger decay heat load of the fuel stored in HSM-15 compared to HSM-1. Considering the 

configurations of HSM–15 and HSM–1 in the module array, the higher temperatures simulated 

for HSM–15 could also reflect the lower exposure of outer surfaces to atmosphere compared to 

HSM–1. However, this effect is expected to be minor compared to the difference in decay heat 

load, which is the principal cause of the higher temperatures realized for HSM–15. The general 

temperature distribution patterns for the HSM–15 and HSM–1 DSC shells and internal surfaces 

are similar. 

Figures 3-20 through 3-23 show temperature distributions for HSM–15 along transverse 

(Figures 3-20 and 3-21) and axial (Figure 3-22 and Figure 3-23) cross sections. The highest 

simulated helium temperature for the HSM–15 DSC internal fuel region in Figure 3-20 {454.48 K 

[359 °F]} is significantly higher than the highest simulated temperature for the HSM–1 DSC 

internal fuel region {390 K [243 °F]}. This is due to the higher decay heat load for HSM–15. 

The spacer plates have the highest temperatures shown in Figure 3-21, as was the case for

HSM–1 (see Figure 3-5), but the maximum temperature of the spacer plates for HSM–15 is 

higher than HSM–1. The concrete at the base remains at a relatively lower temperature 

(almost the same as the ambient temperature) and the pattern is the same for both HSM–1 and 

HSM–15, indicating that the difference in decay heat load did not significantly affect the base 

temperatures. Figures 3-22 and 3-23 demonstrate a sharp temperature difference across the 

DSC shell at the top of the shell.

Axial simulated temperature profiles for the HSM–15 module are presented in Figures 3-24 

through 3-26. While HSM–15 temperatures are generally higher than HSM–1 temperatures, the 

patterns of temperature distribution are similar. Therefore, the HSM–1 observations of

Section 3.2.1.2 apply also to HSM–15.

.
Benchmark Modeling and Sensitivity Analyses

3-15

Figure 3-18. Contours of DSC shell surface temperature for HSM–15 baseline case

Figure 3-19. Temperature distributions for the fuel assembly, spacer plates, and support 

rod for HSM–15 baseline case
Benchmark Modeling and Sensitivity Analyses

3-16

Figure 3-20. Air and helium temperature distributions along transverse cross section for 

HSM–15 baseline case

Figure 3-21. Solid region (concrete, spacer plates) temperature distributions at the 

transverse cross section for HSM–15 baseline case
Benchmark Modeling and Sensitivity Analyses

3-17

Figure 3-22. Temperature distribution for air and helium at axial cross-sectional plane for 

HSM–15 baseline case

Figure 3-23. Temperature distribution for the solid regions at axial cross-sectional plane 

for HSM–15 baseline case
Benchmark Modeling and Sensitivity Analyses

3-18

Figure 3-24. Computed temperature distributions along the axial distance for top and 

right of DSC shell for HSM–15 module

Figure 3-25. Computed temperature distributions along the axial distance for top of heat 

shield and right of heat shield
Benchmark Modeling and Sensitivity Analyses

3-19

Figure 3-26. Computed temperature distributions along the axial distance for fuel region: 

distribution shown for fuel assemblies 1 and 23 for HSM–15 baseline case

3.2.2.2 Velocity Field 

Figures 3-27 through 3-29 present the results for the predicted velocity field of the HSM–15 

module. The airflow path and airflow circulation within the HSM–15 module are similar to those 

for HSM–1. However, the highest magnitude of the velocity {1.18 m/s [3.9 ft/s]} observed in the 

HSM–15 simulations is greater than the corresponding value {1.0 m/s [3.28 ft/s]} in the HSM–1 

simulations observed in Figure 3-7. The higher velocity magnitude can be attributed to the 

higher heat load of HSM–15, which caused stronger convective flow.

3.2.2.3 Comparison with Measured Data

For the HSM–15 module, the only measured temperature is from under the grapple ring. The 

measured temperature is 324 K [124 °F], while the predicted temperature is 335 K [144 °F]. 

Prior to the measurement [Suffield et al 2012], the module door was open for approximately 

20 minutes. For the same reasons discussed for HSM–1 (Section 3.2.1.3), it is expected 

that the measured temperature was lower than the temperature for the normally operated 

closed module.

3.3 Sensitivity Studies

Sensitivity studies of the HSM–1 module simulations were designed to evaluate the effects on 

computed results of mesh refinement, different turbulence models, varying porous media 

resistance, and insolation. In addition, a study was performed to assess the validity of the heat 

transfer coefficient, specified as a boundary condition on the external wall, by performing a 

simulation on an extended domain that included the surrounding atmosphere. Apart from the 
Benchmark Modeling and Sensitivity Analyses

3-20

Figure 3-27. Streamlines showing the airflow path and airflow circulation within the 

HSM–15 module

Figure 3-28. Velocity vector distribution at the transverse cross section for HSM–15 

module
Benchmark Modeling and Sensitivity Analyses

3-21

Figure 3-29. Velocity distribution at the axial cross section for HSM–15 module

sensitivity study on mesh refinement, all simulations used the baseline mesh. The decay 

heat load values calculated for June 2012 were used, as was a 28 °C [301 K, 82.19 °F] 

ambient temperature.

3.3.1 Effect of Mesh Refinement

The objective of this study was to evaluate the sensitivity of calculated temperature values to 

the DSC shell wall y+ parameter, which is defined as the non-dimensional wall distance for a 

wall bounded flow and is often termed the law of the wall coordinate. The y+ parameter is often 

used to describe the quality of mesh in the near wall region. 

The specification of the DSC shell wall y+ parameter controls the location of first node spacing 

perpendicular to the shell wall and the distribution of nodes within the boundary layer. The 

sharp temperature and velocity gradients in the boundary layer provide a strong control on the 

overall heat transfer process and, therefore, temperatures are likely to be especially sensitive to 

inadequate refinement in this zone. The sensitivity study was conducted by adapting and 

refining the mesh near the DSC shell outer wall, but not for the DSC shell inner wall. Due to the 

presence of smaller annuli between the DSC shell and spacer plates, mesh refinement for the 

inner wall was deemed not feasible. It is likely that outer wall mesh refinement has more effect 

on overall heat transfer, since the air flow outside of the DSC shell is turbulent and the helium 

flow inside the DSC shell is laminar. However, for a more comprehensive estimation of the 

effect of y+

, mesh refinement for both inner and outer walls of the DSC shell would be 

necessary. The effect of outer wall mesh refinement on the computed temperature distribution 
Benchmark Modeling and Sensitivity Analyses

3-22

is shown in Figures 3-30 through 3-32. Table 3-2 shows the computational mesh obtained due 

to adaptation and refinement that was used for this sensitivity study and the corresponding 

average y+ value for the first grid point near the DSC outer wall.

Figure 3-30 shows that there is very little difference among the axial distributions of the 

predicted temperatures on the right side of the DSC shell for the three levels of grids. However, 

for the axial temperature distribution on the top of the DSC shell (Figure 3-31), the temperature

distribution for the baseline grid is 5–7 K [9–13 °F] lower than the temperature distributions for 

the two fine discretization levels at the middle of the axial length. This is probably because the 

maximum temperature is predicted in the mid-region at the top of the DSC shell; that region is 

characterized by sharp near-wall gradients. The baseline grid, with a relatively higher value of 

y+

, results in lower values of predicted temperatures near the wall. The baseline grid was not 

fine enough in the near-wall region to capture the flow and temperature gradients in that 

particular region. The two levels of finer grid had enough resolution (lower values of y+

) to 

predict a higher value of the temperature. At the beginning and end of the axial length, the 

temperature distributions for the three grids are almost the same. Figure 3-32 shows an 

insignificant difference in the predicted temperatures in the axial direction for the three grids at 

fuel assembly 1. This may be due to the fact that the mesh was not refined for the DSC internal 

components and, as a result, the heat transfer rates calculated by the three meshes were nearly 

identical. The calculated values of the temperatures for the three levels of refinement are listed 

in Table 3-3. The maximum difference among the results for any location is 1.76 K [3.17 °F] 

and for almost all the observations, the difference is less than 0.5 K [0.9 °F].

In general, the close agreement between the simulated temperatures from all three grid levels 

indicates that all the computational meshes used in the simulations capture the thermal field 

adequately. The majority of the volume had a good quality grid in the original baseline mesh, 

resulting in no significant change in heat transfer and temperature values when the mesh was 

refined. In summary, the simulated temperatures and temperature distributions were generally 

not sensitive to the DSC shell y+

. 

3.3.2 Effect of Turbulence Models

The effect of the choice of turbulence model on the computed temperature distribution is shown 

in Figures 3-33 through 3-36. The turbulence models considered in the sensitivity study were:

(i) standard k-ε model with scalable wall function, 

(ii) k-ε model with near-wall low Reynolds number correction and buoyancy effects without 

wall functions (baseline case), 

(iii) standard shear stress transport (SST) k-ω model with low Reynolds number correction 

without wall functions, and

(iv) SST k-ω model with buoyancy and low Reynolds number correction without 

wall functions. 

Figures 3-33 and 3-34 show the axial temperature distributions on the right side of the DSC 

shell and on the top of the DSC shell for the turbulence models. Figures 3-35 and 3-36 show 

the axial temperature distribution in fuel assemblies 1 and 23, respectively, for the turbulence 

models. The simulated results for the SST k-ω models (both with and without buoyancy) are 

almost identical. For the k-ε model, results show that including low Reynolds number correction 

with buoyancy resulted in lower values of the predicted temperatures than any other model, and 
Benchmark Modeling and Sensitivity Analyses

3-23

Figure 3-30. Computed temperature distributions along the axial distance for right side 

of DSC shell: mesh refinement for HSM–1 module

Figure 3-31. Computed temperature distributions along the axial distance for top of DSC 

shell: mesh refinement for HSM–1 module
Benchmark Modeling and Sensitivity Analyses

3-24

Figure 3-32. Computed temperature distributions along the axial distance for fuel 

assembly 1: mesh refinement for HSM–1 module

Table 3-2. Computational mesh used in the mesh refinement study

Simulation 

cases for 

HSM–1

Total 

number of 

grids

Minimum y+ on 

DSC outer wall

Maximum y+ on 

DSC outer wall

Surface Average 

y+ on DSC outer 

wall

Baseline Grid 3.2 × 106 0.054 12.21 2.97

Refinement-1 4.2 × 106 0.038 10.60 1.575

Refinement-2 6.5 × 106 0.027 7.01 0.8193

Table 3-3. Comparison of predicted temperatures for the three different levels of mesh: 

simulations for HSM–1 module

Location

Baseline mesh (K)

[°F = 1.8 K-460]

Fine mesh 1

(refinement-1) (K)

[°F = 1.8 K-460]

Fine mesh 2

(refinement-2) (K)

[°F = 1.8 K-460]

Under Grapple ring 324.18 325.29 324.20

0-m-side-right 323.04 323.70 322.94

0-m-side-top 325.18 326.04 324.98

0-m-side-left 322.68 323.61 322.81

0-m-rail-right 319.05 319.09 319.51

0-m-rail-left 319.05 320.01 319.67

0.51-m-side-right 325.74 326.94 327.17

0.51-m-side-top 336.94 336.01 336.25

0.51-m-side-left 324.94 326.04 326.63

0.51-m-rail-right 318.81 319.31 319.48
Benchmark Modeling and Sensitivity Analyses

3-25

Table 3-3. Comparison of predicted temperatures for the three different levels of mesh: 

simulations for HSM–1 module

Location

Baseline mesh (K)

[°F = 1.8 K-460]

Fine mesh 1

(refinement-1) (K)

[°F = 1.8 K-460]

Fine mesh 2

(refinement-2) (K)

[°F = 1.8 K-460]

0.51-m-rail-left 318.7 319.16 319.29

1.02-m-side-right 334.35 334.05 334.51

1.02-m-side-top 346.89 348.65 347.68

1.02-m-side-left 333.39 333.86 334.34

1.02-m-rail-right 322.04 322.22 322.74

1.02-m-rail-left 321.97 321.93 322.44

Figure 3-33. Computed temperature distributions along the axial distance for right side 

of DSC shell: effect of turbulence model for HSM–1 module
Benchmark Modeling and Sensitivity Analyses

3-26

Figure 3-34. Computed temperature distributions along the axial distance for top of DSC 

shell: effect of turbulence models for HSM–1 module

Figure 3-35. Computed temperature distributions along the axial distance for fuel 

assembly 1: effect of turbulence model for the HSM–1 module
Benchmark Modeling and Sensitivity Analyses

3-27

Figure 3-36. Computed temperature distributions along the axial distance for fuel 

assembly 23: effect of turbulence model for the HSM–1 module

the temperatures are closest to the measured values at all locations. The calculated values of 

the temperatures for all the turbulence models are listed in Table 3-4. In general, the calculated 

values are numerically close, with the maximum difference among predicted temperatures being 

10.15 K [18.27 °F], at the top central part of the DSC shell.

3.3.3 Effect of Porous Media Resistance 

The effects of porous media flow resistance on the computed temperature distribution are 

shown in Figures 3-37 through 3-40. The cases considered include the baseline case, a higher 

porous media resistance compared to the baseline case, a lower porous media resistance 

compared to the baseline case, zero resistance, and one case with no porous media. In the 

zero resistance case, the porous media settings of the solver were retained, but the individual 

resistance coefficients were set to zero. For the no porous media case, the porous media 

settings within the solver were turned off and the fuel assembly region was treated as a regular 

fluid zone. The higher and lower porous media resistance coefficients were arbitrarily set as 

one order of magnitude higher and lower than the baseline case. Table 3-5 lists the viscous 

and resistance coefficients that were used in this study.

Figures 3-37 and 3-38 show the axial temperature distribution on the right side of the DSC 

shell and on the top of DSC shell for the different values of the porous media resistance. 

Figures 3-39 and 3-40 show the axial temperature distribution for fuel assemblies 1 and 23, 

respectively, for the different porous media resistance values. The simulation results show that 

increasing or decreasing the porous media resistance relative to the baseline value does not 

significantly affect the temperature distribution. 
Benchmark Modeling and Sensitivity Analyses

3-28

Table 3-4. Comparison of predicted temperatures for different turbulence models: 

simulations for HSM–1 module

Location

Baseline 

case, low 

re standard 

k-ε with 

buoyancy

(K)

[°F = 1.8 K460]

Standard k-ε 

model with 

scalable wall 

function and 

buoyancy 

effects (K)

[°F = 1.8 K460]

SST k-ω model 

with low 

Reynolds 

number 

corrections (K)

[°F = 1.8 K-460]

SST k-ω model 

with low Reynolds 

number corrections 

and buoyancy 

effects (K)

[°F = 1.8 K-460]

Under grapple ring 324.18 327.99 327.02 326.922

0-m-side-right 323.04 326.92 325.63 325.64

0-m-side-top 325.18 328.82 328.65 328.39

0-m-side-left 322.68 326.78 325.35 321.21

0-m-rail-right 319.05 322.75 321.15 321.26

0-m-rail-left 319.05 322.56 321.05 320.98

0.51-m-side-right 325.74 332.17 329.46 329.83

0.51-m-side-top 336.94 342.38 342.64 342.51

0.51-m-side-left 324.94 332.04 328.83 329.04

0.51-m-rail-right 318.81 323.27 320.52 320.64

0.51-m-rail-left 318.7 323.29 320.21 320.64

1.02-m-side-right 334.35 342.41 338.3 338.68

1.02-m-side-top 346.89 356.68 357.04 355.74

1.02-m-side-left 333.39 342.19 337.29 337.33

1.02-m-rail-right 322.04 327.22 323.25 323.38

1.02-m-rail-left 321.97 327.2 322.9 323.65

Figure 3-37. Computed temperature distributions along the axial distance for right side 

of DSC shell: effect of porous media for HSM–1 module
Benchmark Modeling and Sensitivity Analyses

3-29

Figure 3-38. Computed temperature distributions along the axial distance for top of DSC 

shell: effect of porous media for HSM–1 module

Figure 3-39. Computed temperature distributions along the axial distance for fuel 

assembly 1: effect of porous media for HSM–1 module
Benchmark Modeling and Sensitivity Analyses

3-30

Figure 3-40. Computed temperature distributions along the axial distance for fuel 

assembly 23: effect of porous media for HSM–1 module

Table 3-5. Resistance coefficients used for porous media resistance 

parametric study

Baseline

High 

resistance

Low 

resistance

Zero 

resistance

No porous 

media

Viscous 

Resistance 

(1/m2 )

[1/m2 =

0.0931/ft2

]

Axial 

Direction

1.284 ×

105

107 1000 0.0 Not 

specified

Radial 

Direction-1

108 108 108 0.0 Not 

specified

Radial 

Direction-2

108 108 108 0.0 Not 

specified

Inertial 

Resistance 

(1/m) 

[1/m = 0.305 

1/ft]

Axial 

Direction

0.78 10 0.01 0.0 Not 

specified

Radial 

Direction-1

100 100 100 0.0 Not 

specified

Radial 

Direction-2

100 100 100 0.0 Not 

specified

However, the zero resistance and no porous media cases do alter the basic temperature 

distribution pattern, particularly for fuel assemblies 1 and 23. The calculated values of the 

temperatures and their comparisons for all the different porous media resistance cases are 

shown in Table 3-6. The values listed in the table substantiate the observation that there is no 

significant difference among the baseline, high resistance, and low resistance cases, with the 
Benchmark Modeling and Sensitivity Analyses

3-31

Table 3-6. Comparison of predicted temperatures (K) for different porous media 

resistance values: simulations for HSM–1 module

Location

Baseline

(K)

[°F =

1.8 K-460]

High 

resistance

(K)

[°F =

1.8 K-460]

Low 

resistance

(K)

[°F =

1.8 K-460]

Zero 

resistance

(K)

[°F =

1.8 K-460]

No porous 

media (K)

[°F =

1.8 K-460]

Under grapple 

ring

324.18 324.88 324.77 333.12 328.33

0-m-side-right 323.04 323.57 323.48 330.94 327.09

0-m-side-top 325.18 325.68 325.58 334.93 330.08

0-m-side-left 322.68 323.27 323.19 330.40 326.68

0-m-rail-right 319.05 319.88 319.76 324.13 321.75

0-m-rail-left 319.05 319.69 319.56 323.86 321.75

0.51-m-side-right 325.74 326.05 325.85 327.47 328.56

0.51-m-side-top 336.94 335.79 335.46 337.68 339.35

0.51-m-side-left 324.94 325.38 325.34 327.45 327.92

0.51-m-rail-right 318.81 319.39 319.23 319.92 319.61

0.51-m-rail-left 318.7 319.14 318.94 319.68 319.35

1.02-m-side-right 334.35 334.47 334.26 331.40 332.07

1.02-m-side-top 346.89 347.31 347.51 342.85 343.92

1.02-m-side-left 333.39 333.95 333.82 330.81 331.74

1.02-m-rail-right 322.04 322.53 322.34 320.39 320.44

1.02-m-rail-left 321.97 322.25 322.01 320.04 320.04

maximum difference among predicted results being 1.5K [2.7 °F]. But the zero and no porous 

media approximations yield a significantly different temperature pattern. It is expected that the 

zero resistance and no porous media specification cases would yield identical outcomes. 

However, the difference in the patterns of temperature distributions between these two cases 

indicates that the source terms in the flow solver specified for porous media with zero resistance 

may have affected the solution, though these sources are expected to be zero. In general, this 

study shows that the calculated temperatures are not sensitive to the exact resistance 

coefficients, but the fuel assemblies cannot be assumed as free flow zones with no 

porous resistance.

3.3.4 Effect of Insolation

The effect of insolation was evaluated using two cases: the baseline case (which includes 

insolation) and the baseline case modified to omit insolation. The solar insolation loads in the 

baseline case (see Section 2.4) are taken from the values prescribed in the regulations for 

packaging and transportation of radioactive material [10 CFR Part 71]. Section 71.71 of 

10 CFR Part 71 specifies a solar load of800 g-cal/cm2 [2,949.35 BTU/ft2

] for horizontal surfaces 

and 200 g-cal/cm2 [737.337 BTU/ft2

] for flat non-horizontal surfaces applied over a 12-hour 

period. For the simulation case without insolation, the simulation settings and the parameters 

are identical to the baseline case except that no insolation flux is specified. Figure 3-41 shows 

the axial temperature distribution on the top of DSC shell and on the right side of the DSC shell 

for the two cases. Figure 3-42 shows the axial temperature distribution for fuel assemblies 1 

and 23 for the two cases. Both Figures 3-41 and 3-42 indicate that omitting insolation reduces 

the simulated temperatures, which is expected because insolation provides an additional heat 
Benchmark Modeling and Sensitivity Analyses

3-32

Figure 3-41. Computed temperature distributions along the axial distance for right side 

and top of DSC shell: effect of insolation for HSM–1 module

Figure 3-42. Computed temperature distributions along the axial distance for fuel 

assemblies 1 and 23: effect of insolation for HSM–1 module
Benchmark Modeling and Sensitivity Analyses

3-33

load on the external surfaces of the module. Peak temperatures are affected only by 

approximately 1 to 2 K [2 to 4 °F] by adding insolation, but the temperature patterns differ 

slightly and temperatures nearer the DSC ends differ by up to 4 K [7.2 °F]. 

The calculated values of the temperatures for the two cases are shown at the measurement 

points from Table 3-1 in Table 3-7. The values in the table substantiate that including insolation 

consistently increases predicted temperatures. The maximum temperature difference between 

the results obtained using insolation and the baseline case is 3.77 K [6.78 °F]. In general, this 

study shows that insolation does affect the internal temperature of the storage system and 

therefore should be considered.

3.3.5 Effect of Including the Surrounding Atmosphere in the 

Computational Domain

A sensitivity study was performed to analyze the effect of an extended domain on the predicted 

temperature field. The primary objective of this specific study was to assess the validity of the 

heat transfer coefficient specified on the outer walls for the baseline cases, which use a 

computational domain that does not include the extended atmospheric domain. The baseline 

case uses a constant heat transfer coefficient of 10 W/m2

K [1.70 BTU h−1 ft−2 °F−1] applied on 

external surfaces; in the sensitivity study, surface-averaged heat transfer coefficients were 

derived from simulation results. For this study with extended domain, the standard k-ε model 

with scalable wall function was used. This choice of turbulence model was based on the quality 

of the grid on the outer surface of the concrete module, which was not as fine as the mesh 

around the DSC shell. At the outer exterior of the domain, a pressure inlet boundary condition 

with zero gauge pressure was specified. The bottom surface of the domain was specified as a 

wall and the back surface of the domain was specified to have a symmetric boundary condition.

Figure 3-43 shows the extended domain used for this study along an axial plane. The domain 

extends to 2.5 times the height of the HSM–1 module in the vertical direction, 1.5 times the 

length of the HSM–1 module in the horizontal direction, and 1.5 times the width of the HSM–1 

module in the span wise directions. The boundaries of the extended domain are maintained at 

ambient temperature and pressure. The sensitivity study is expected to slightly underestimate 

heat transfer coefficients because wind, which acts to cool the module, was not included.

Figures 3-44 and 3-45 show simulated temperatures and velocity magnitudes at an axial 

mid-plane of the computational domain. Temperatures calculated from the extended domain 

simulation were compared with computations performed on the regular domain using the same 

turbulence model (Section 3.3.2). The temperature differences between these two cases in the 

fuel assembly region were less than 2.5 K [4.5 °F]. The velocity distribution and magnitude 

within the module are similar to those of the baseline case observed in Figure 3-13; outside the 

module, the heated air rises vertically. Figure 3-46 shows the computed convective heat 

transfer coefficient distribution on the external concrete surface. Using simulation results with 

the extended atmospheric domain to calculate surface-averaged heat transfer coefficients 

results in values of 9.467 W/m2

K [1.67 BTU h−1 ft−2 °F−1] and 9.240 W/m2K [1.63 BTU h−1 ft−2 

°F−1] for horizontal and vertical surfaces, respectively. These compare well with the value of 

10 W/m2

K [1.70 BTU h−1 ft−2 °F−1] used in the baseline simulations and would be expected to 

increase slightly if wind were included in the study. Accordingly, the heat transfer coefficient 

used for the baseline case appears to be a reasonable value.
Benchmark Modeling and Sensitivity Analyses

3-34

Table 3-7. Comparison of predicted temperatures for with and without insolation: 

simulations for HSM–1 module

Location

Baseline (K)

[°F = 1.8 K-460]

Without insolation (K)

[°F = 1.8 K-460]

Under Grapple ring 324.18 320.23

0-m-side-right 323.04 320.11

0-m-side-top 325.18 321.95

0-m-side-left 322.68 319.83

0-m-rail-right 319.05 316.32

0-m-rail-left 319.05 316.28

0.51-m-side-right 325.74 323.92

0.51-m-side-top 336.94 334.09

0.51-m-side-left 324.94 323.19

0.51-m-rail-right 318.81 316.70

0.51-m-rail-left 318.70 316.75

1.02-m-side-right 334.35 332.96

1.02-m-side-top 346.89 344.85

1.02-m-side-left 333.39 331.75

1.02-m-rail-right 322.04 320.23

1.02-m-rail-left 321.97 320.41

Figure 3-43. Extended domain for the study of the effect of atmosphere for 

HSM–1 module
Benchmark Modeling and Sensitivity Analyses

3-35

Figure 3-44. Temperature contours at an axial mid-plane of the extended domain: effect 

of atmosphere for HSM–1 module

Figure 3-45. Velocity contours at an axial mid-plane of the extended domain: effect of 

atmosphere for HSM–1 module
Benchmark Modeling and Sensitivity Analyses

3-36

Figure 3-46. Computed heat transfer coefficient distribution over the HSM–1 module: 

effect of extended domain and atmosphere

3.4 Observations From Benchmark Study

Detailed 3-D computational models were developed for thermal analyses of the HSM–1 and 

HSM–15 modules. Simulations were performed to benchmark the developed models by 

comparing the numerical predictions with measured data and to investigate the sensitivity of the 

results to model assumptions and boundary conditions. The major observations from this 

benchmark study are as follows:

• The simulations provided valuable insight into the heat transfer process inside the DSC. 

The buoyancy-driven natural convection flow of the air and the helium gas influences the 

fuel region temperature and the DSC shell temperature. 

• Predicted temperatures were consistently higher than values measured in the field. The 

measured values are subject to significant uncertainty due to a number of factors,

including difficulty in probe insertion, open module doors that permitted atmospheric 

mixing, and a one-time measurement without any repetition. Therefore, the measured 

temperatures may not be truly representative of typical operating conditions with a 

closed door.

• A parametric study of sensitivity of the computed results to y+ and mesh refinement to 

improve wall y+ showed no significant change in overall heat transfer and temperature 

distribution when the mesh is refined relative to the baseline mesh.

• The baseline low Reynolds k-ε model-predicted temperatures were closest to the 

measured values. The maximum and minimum deviations for this case were 25.7 K 

[46.26 °F] and 3.33 K [5.99 °F], respectively.
Benchmark Modeling and Sensitivity Analyses

3-37

• A sensitivity study with the turbulence model showed that including the low Reynolds 

number correction in the k-ε model resulted in lower values of the predicted 

temperatures on the DSC outer shell.

• An order of magnitude increase or decrease in porous media resistance (compared to 

baseline) did not affect the temperature distribution. However, zero resistance and 

assuming no porous media significantly flattened the temperature distribution pattern.

• Omitting the solar flux on external surfaces resulted in a temperature drop of up to 4 K 

[7.2 °F] and steeper temperature gradients, although peak temperatures were affected 

by a smaller amount.

• The predicted temperatures with an extended domain were comparable to those

obtained from the regular domain simulation using the same turbulence model. The 

surface-averaged heat transfer coefficient values used in the baseline case were found 

to be reasonable based on the extended domain study.

4-1

4. THERMAL EVOLUTION SIMULATIONS

This chapter presents the results of simulations performed to predict storage module 

temperatures for up to 300 years of storage in horizontal storage modules HSM–1 and 

HSM–15. The analysis was performed using a sequence of steady-state simulations based on 

the projected decay heat over the 300-year period. 

4.1 Numerical Approach to Thermal Evolution Analysis 

A multi-step, steady-state analysis approach was used to simulate the thermal behavior of the 

storage modules. In this approach, a steady-state analysis was performed at selected time 

points covering the 300-year time span. At each time step, quasi-steady behavior of the flow 

and thermal field within and around the storage modules was assumed. Independently 

calculated heat loads at the discrete time points were used in the simulations.

For both HSM–1 and HSM–15, a 5-year interval was used to calculate the discrete time points 

for the first 100 years after storage, while a 10-year interval was used for the remaining 

200 years of storage. Simulations were specifically performed for the calendar year 2012 

because the observations were made at that time for both HSM-1 and HSM-15. For HSM–1, 

the fuel was loaded in the storage module in 1993. To span the 300 years of storage, 

simulations in this study were performed at specified times from 1993 to 2293. For the HSM–15 

module, the fuel was loaded in year 1996 and simulations were conducted until year 2297, to 

cover the 301-year storage time. The span of simulation was extended by a year for HSM-15 to 

accommodate the calendar year 2012 in the series and maintain time steps of 5 and 10 years. 

The heat load values specified at discrete time points were based on the decay characteristics 

of the fuel and inventory, calculated in a separate analysis. The calculated decay heat load 

values used in the simulations for some specific years are provided in Tables 4-1 and 4-2 for

HSM–1 and HSM–15, respectively. The fuel assembly loading identification and locations are 

shown in Figures 2-13 and 2-14 for HSM–1 and HSM–15, respectively. The heat load 

calculated for each of the fuel assemblies at different years was applied as a uniform volumetric 

heat generation rate over the active fuel length within each fuel region. All simulations 

performed for this analysis used the baseline computational grid and input parameters for both 

HSM–1 and HSM–15 (see Chapter 2 for a description). Insolation was included and the k-ε 

turbulence model with low Reynolds number correction and buoyancy effects was used. 

Selected simulation results are presented for calendar years 1993, 2093, 2193, and 2293 for

HSM–1, and calendar years 1996, 2097, 21967 and 2297 for HSM–15, to highlight the evolution 

of temperature with time. Chapter 3 presents results for the year 2012 for both 

HSM–1 and HSM–15 as the baseline simulations. 

4.2 Results of Thermal Evolution Simulations 

The thermal evolution simulations focus on the maximum and minimum temperatures for 

particular module components, including the cladding, fuel sleeves, spacer plates, the dry 

storage canister (DSC) shell, and the heat shield. Snapshots along a transverse cross section 

for the temperature distribution in the airflow pathways and within the DSC are also provided for 

context. The maximum and minimum temperatures in the fuel assembly region are used to 

represent maximum and minimum cladding temperatures, because the fuel assembly region is 

approximated as a continuous porous region (see Chapter 2).
Thermal Evolution Simulations

 4-2

Table 4-1. Calculated decay heat load values at different years for HSM–1 module

Assembly 

ID

Calculated 

1993 heat 

load

(start of 

storage)

(Watt)

1 Watt =

[3.413 

BTU/hr]

2012 heat load

(current year of 

measurement)

(Watt)

1 Watt =

[3.413 BTU/hr]

Calculated 

2093 heat 

load

(100 years of 

storage)

(Watt)

1 Watt =

[3.413 

BTU/hr]

Calculated 

2193 heat 

load

(200 years of 

storage)

(Watt)

1 Watt =

[3.413 

BTU/hr]

Calculated 

2293 heat 

load

(300 years of 

storage)

(Watt)

1 Watt =

[3.413 

BTU/hr]

1A007 181.33 130.5 50.67 33 27.767

1A034 181.7 130.7 50.77 33 27.233

1A035 184.67 133.0 51.7 33.7 28.367

1A023 190.63 137.5 53.6 34.9 29.367

1A041 190.47 137.4 53.5 34.9 29.367

1A059 190.43 137.4 53.5 34.9 29.367

1A029 191.4 138.2 53.8 35.1 28.933

1A054 192.9 139.3 54.3 35.4 29.767

1B032 194.9 139.9 52.8 33.7 28.267

1B013 266.3 183.9 71.3 44.9 36.767

1B053 267.37 184.7 71.5 45.067 36.867

1B079 271.37 187.4 72.6 45.7 37.367

1B046 278.37 192.2 74.4 46.767 38.267

1B044 283.5 195.5 75.7 47.567 38.867

1B004 286.67 197.8 76.5 48.067 39.267

1B018 286.67 197.7 76.5 48 39.267

1B061 287 198 76.6 48.067 39.267

1B003 288.23 198.9 76.9 48.267 39.467

1C005 254.73 175.2 65.1 40.3 32.967

1C103 275.76 189.9 71.3 44.1 35.967

1C107 275.07 189.4 71.1 44 35.867

1C106 287.8 198.3 74.5 46.067 37.567

1C206 296.6 202.1 76.8 47.467 38.67

1C214 296.1 201.7 76.7 47.367 38.567
Thermal Evolution Simulations

 4-3

Table 4-2. Calculated decay heat load values at different years for HSM–15 module

Assembly 

ID

Calculated 

1996 heat 

load

(start of 

storage)

(Watt)

1 Watt =

[3.413 

BTU/hr]

2012 heat load

(current year of 

measurement)

(Watt)

1 Watt =

[3.413 BTU/hr]

Calculated 

2097 heat 

load

(101 years of 

storage)

(Watt)

1 Watt =

[3.413 

BTU/hr]

Calculated 

2197 heat 

load

(201 years of 

storage)

(Watt)

1 Watt =

[3.413 

BTU/hr]

Calculated 

2297 heat 

load

(301 years of 

storage)

(Watt)

1 Watt =

[3.413 

BTU/hr]

2F177 346.1 272 83.83 50.46 40.75

2F136 367.3 291 89.09 53.52 43.04

1E 101 404.60 301 98.27 59.79 48.07

2F021 533.97 385 120.90 70.60 55.63

2F024 533.90 385 120.90 70.59 55.63

2F118 360.3 286 87.40 52.53 42.30

2F179 367.1 290 89.04 53.49 43.02

1H115 545.20 406 119.97 69.46 54.48

2G123 568.80 427 123.93 71.43 55.83

1H109 545.80 407 120.07 69.50 54.55

1G018 445.33 327 105.10 62.41 49.74

2F117 359.0 285 87.06 52.33 42.10

2F123 359.3 285 87.15 52.39 42.19

1G027 456.97 340 107.57 63.74 50.68

1H119 547.60 408 120.40 69.68 54.68

1H124 546.10 407 120.10 69.53 54.57

1E 102 404.50 301 98.24 59.77 48.06

1F119 302.93 232 78.96 49.44 40.34

2F166 347.9 272 84.92 51.28 41.38

1E 121 406.60 302 98.70 60.04 48.27

1E 120 406.27 302 98.62 59.99 48.23

1F110 291.93 225 76.11 47.71 39.05

2F125 345.4 274 83.70 50.37 40.65

1F113 294.97 225 77.45 48.67 39.86
Thermal Evolution Simulations

4-4

The minimum and maximum temperatures decline over the 300 years of storage for all 

components (Figures 4-1 through 4-5). All components show a consistent trend, with a 

relatively rapid rate of temperature drop over the first 100 years and a more gradual rate of 

change thereafter. For example, the maximum cladding temperature decreases by 

approximately 129 K [232 °F] within the first 100 years and 28 K [50.4 °F] during the remaining 

200 years of storage. This decline is driven by the change in decay heat load. The maximum 

component temperatures show much larger changes than the minimum component 

temperatures, which are relatively close to the ambient temperature throughout the storage 

period for all components except the cladding. Most of the minimum component temperatures 

decline by 10 K [18 °F] or less, reaching an asymptotic value only a few degrees above the 

ambient temperature partway through the storage period. Maximum component temperatures 

consistently drop with distance from the fuel rods (i.e., cladding to fuel sleeve to spacer plates to 

DSC shell to heat shield), as do the minimum temperatures, but at all times the temperature 

difference between components is much larger for the maximum temperature than the 

minimum temperature.

Figures 4-6 through 4-9 show the evolution of the temperature distribution in the airflow 

pathways and within the DSC along a transverse plane in HSM–1 for the same calendar years

shown in Figures 4-1 through 4-5. The relative spatial temperature distribution remains 

consistent from century to century, but the most rapid decline occurs in the high 

temperature zone in the DSC. The peak temperature in the fuel region decreases from 370 to 

320 K [207 to 116 °F] over the 300 years of storage, consistent with Figures 4-1 and 4-4, and 

the temperature change is greater in the first 100 years than the last 200 years of storage.

Figures 4-10 through 4-18 present the same simulated results for HSM–15 that are shown in 

Figures 4-1 through 4-9 for HSM–1. There is a similar pattern of temperature change in HSM–1 

and HSM–15, except that the simulated temperatures for the HSM–15 module are higher than 

for HSM–1. Higher temperatures are expected in HSM–15 mainly because the decay heat load 

is nearly twice as large in HSM–15 (compare Tables 4-1 and 4-2). The corner position of 

HSM–1 in the module pattern exposes a side to atmospheric cooling, while HSM–15 is 

surrounded by other modules, but the effect of location is minor in comparison to the effect of 

the higher decay heat load.

Tables 4-3 and 4-4 present the calculated maximum and minimum cladding and DSC shell 

temperatures at periodic intervals for the HSM–1 and HSM–15 modules, respectively. 

These tables confirm that temperatures change most rapidly over the first 100 years of 

storage and then asymptotically converge to values only mildly above the ambient temperature 

of 301 K [82 °F].
Thermal Evolution Simulations

4-5

Figure 4-1. Variation of maximum and minimum cladding temperature over 300 years for 

HSM–1 module

Figure 4-2. Variation of maximum and minimum fuel sleeve temperature over 300 years 

for HSM–1 module 
Thermal Evolution Simulations

4-6

Figure 4-3. Variation of maximum and minimum spacer plates temperature over 300 

years for HSM–1 module

Figure 4-4. Variation of maximum and minimum DSC shell temperature over 300 years 

for HSM–1 module
Thermal Evolution Simulations

4-7

Figure 4-5. Variation of maximum and minimum heat shield temperature over 300 years 

for HSM–1 module

Figure 4-6. Air and helium temperature distributions at a transverse plane for HSM–1:

year 1993
Thermal Evolution Simulations

4-8

Figure 4-7. Air and helium temperature distributions at a transverse plane for HSM–1: 

year 2093 (100 years after storage)

Figure 4-8. Air and helium temperature distributions at a transverse plane for HSM–1: 

year 2193 (200 years after storage)
Thermal Evolution Simulations

4-9

Figure 4-9. Air and helium temperature distributions at a transverse plane for HSM–1: 

year 2193 (300 years after storage)

Figure 4-10. Variation of maximum and minimum cladding temperature over 301 years 

for HSM–15 module
Thermal Evolution Simulations

4-10

Figure 4-11. Variation of maximum and minimum fuel sleeve temperature over 301 years 

for HSM–15 module

Figure 4-12. Variation of maximum and minimum spacer plate temperature over 301

years for HSM–15 module
Thermal Evolution Simulations

4-11

Figure 4-13. Variation of maximum and minimum DSC shell temperature over 301 years

for HSM–15 module 

Figure 4-14. Variation of maximum and minimum heat shield temperature over 301 years 

for HSM–15 module
Thermal Evolution Simulations

4-12

Figure 4-15. Air and helium temperature distributions at a transverse plane for HSM–15: 

year 1996

Figure 4-16. Air and helium temperature distributions at a transverse plane for HSM–15: 

year 2097 (101 years after storage)
Thermal Evolution Simulations

4-13

Figure 4-17. Air and helium temperature distributions at a transverse plane for HSM–15: 

year 2197 (201 years after storage)

Figure 4-18. Air and helium temperature distributions at a transverse plane for HSM–15: 

year 2297 (301 years after storage)
Thermal Evolution Simulations

4-14

Table 4-3. Calculated temperature values at different years for HSM–1 module

Year

Interval 

(years)

Maximum 

cladding 

temperature (K)

°F = [1.8K-460]

Minimum 

cladding 

temperature (K)

°F = [1.8 K-460]

Maximum DSC 

shell 

temperature (K)

°F = [1.8 K-460]

Minimum DSC 

shell 

temperature (K)

°F = [1.8 K-460]

1993 0 414.392 317.645 371.69 311.62

2012 19 386.316 313.742 354.01 308.81

2032 39 369.935 311.472 345.19 307.47

2052 59 357.009 309.686 337.28 306.29

2072 79 348.25 308.61 331.54 305.66

2093 100 341.366 307.719 327.55 305.09

2123 130 335.907 306.813 324.29 304.42

2153 160 331.841 306.121 322.51 304.36

2183 190 330.095 305.991 320.76 303.82

2213 220 328.467 305.812 319.67 303.77

2243 250 327.206 305.59 318.95 303.64

2273 280 326.179 305.489 318.57 303.57

2293 300 324.832 305.077 317.91 303.47
Thermal Evolution Simulations

 4-15

Table 4-4. Calculated temperature values at different years for HSM–15 module

Year

Interval 

(years)

Maximum 

cladding 

temperature (K)

Minimum 

cladding 

temperature (K)

Maximum DSC 

shell temperature 

(K)

Minimum DSC 

shell temperature 

(K)

1996 0 490.834 324.149 412.74 314.83

2017 21 438.218 318.443 381.41 311.27

2037 41 407.502 314.087 363.79 308.99

2057 61 386.23 312.261 351.66 307.45

2077 81 371.55 310.47 343.44 306.37

2097 101 361.442 309.194 337.65 305.61

2127 131 351.688 307.962 332.06 304.85

2157 161 345.589 307.183 328.74 304.35

2187 191 341.575 306.676 326.56 303.99

2217 221 338.803 306.288 324.86 303.79

2247 251 336.706 305.902 323.54 303.57

2277 281 335.008 305.66 322.43 303.42

2297 301 333.928 305.561 322.01 303.35
Thermal Evolution Simulations

4-16

4.3 Observations from the Thermal Evolution Study

Detailed three-dimensional numerical simulations were performed to predict the thermal 

behavior of the HSM–1 and HSM–15 modules for 300 years of storage time. The analysis used 

steady-state simulations at multiple discrete time points, each corresponding to late June in a 

selected calendar year. The major observations from this study are the following:

• The simulations provided valuable insight into the time-dependent variation of maximum 

and minimum temperatures for cladding and internal components such as fuel sleeves, 

spacer plates, the DSC shell, and the heat shield. The general pattern of temperature 

change is similar for all the components, but the degree of variation is different for 

different components.

• Simulated results for both HSM–1 and HSM–15 show that the maximum temperature of 

any component drops substantially within the first 100 years of storage, with subsequent 

changes in temperature being more gradual.

• Simulated results highlight that the temporal variation of minimum component 

temperatures is small, approximately 10 K [18 °F].

• The general patterns of the temperature changes for HSM–1 and HSM–15 are similar. 

As with the baseline cases considered in Section 3, simulated temperatures for HSM–15 

are uniformly higher than HSM–1 because of the larger decay heat load of the fuel 

contained in HSM–15 and the greater exposure of HSM–1 to the atmosphere.
5-1

5. UNCERTAINTY QUANTIFICATION AND GRID CONVERGENCE 

INDEX (GCI) ANALYSIS

The credibility and predictive capability of computational fluid dynamics (CFD) is greatly 

enhanced, and CFD can be considered a reliable tool for engineering analysis, if the numerical 

errors and uncertainties present in the solution can be reliably estimated [Phillips and Roy 2011;

Rezende et al. 2012; Oberkampf and Roy 2010]. Error and uncertainty for a solver and a 

solution is formally estimated through a verification and validation (V&V) process. In 2009, the 

American Society of Mechanical Engineers (ASME) published a standard [ASME 2009]

explaining the procedures for code verification, solution verification, and validation of CFD 

studies. The general objective of the V&V process is to estimate the accuracy of a 

numerical study. 

Validation is the process used to determine whether a model accurately represents the physical 

phenomena of interest and to obtain an estimate of modeling uncertainty. This is primarily done 

by comparing simulated results with suitable experimental data. Verification deals with the 

uncertainty of the solution that can be attributed to numerical error, irrespective of how well the 

solution represents the actual physical process it models. A verification study consists of two 

components: code verification and solution verification. Code verification is used to 

demonstrate that the mathematical equations are correctly programmed. It ensures that the 

solver is free of any error due to coding and that the solution algorithm has been implemented 

correctly. In general, developers perform code verification based on evaluation of a known

problem by the solver. In contrast, solution verification does not deal with the correctness of a 

code, but involves estimation of the numerical error or uncertainty for specific simulation results

and is performed by users of a code. Discretization errors and iteration errors are the two main 

sources of numerical error. The solution verification process is also known as numerical error 

estimation [Celik et al. 2008]. 

This chapter presents the methodology and results from a numerical error estimation analysis 

performed for the dry cask thermal analysis described in previous chapters. The solution 

verification process described in ASME V&V 20-2009 [ASME 2009] was used to quantify 

numerical uncertainty for the calculated temperatures in horizontal storage modules HSM–1 and 

HSM–15. The uncertainty analysis was performed for solution verification with systematic grid 

refinement. For this study, two representative variables—maximum and minimum cladding 

temperature—were used to indicate solution behavior in addressing solution verification. These 

variables provide convergence criteria for numerical error estimation that are directly applicable 

to understanding the performance of the entire storage system. Although in principle the spatial 

location of these representative variables might shift from grid to grid, in practice the locations 

were observed to be essentially invariant. Therefore, the maximum and minimum cladding 

temperatures were used as representative parameters for numerical error estimation.

5.1 Uncertainty Analysis: Definitions and Concepts

The V&V process and concepts are briefly described in this section to facilitate the discussion of 

results, which follow.

5.1.1 Code Validation

The objective of code validation is to estimate (i) how accurately a model represents the physics 

of a problem and (ii) estimate the modeling error within a specified uncertainty range 

[ASME 2009]. Generally, this is achieved by comparing the result of a given simulation with an 
Uncertainty Quantification and Grid Convergence Index (GCI) Analysis

5-2

experiment. Because code validation involves both numerical modeling results and 

experimental data, the source of uncertainty can be due to model inputs, numerical issues 

related to the model, or experimental uncertainty. As shown in ASME V&V 20-2009

[ASME 2009], if all these sources of errors are independent, then validation uncertainty can be 

defined as

Equation (5-1)

where

uval = Validation uncertainty

unum = Numerical uncertainty

uinput = Input uncertainty

uexp = Experimental uncertainty.

Numerical uncertainty contributes to the overall assessment of validation and is not the entire 

measure of the total uncertainty. Numerical uncertainty is estimated through verification

(see Section 5.1.2). The present study considers numerical uncertainty without considering 

input or experimental uncertainty. To accurately measure the total uncertainty associated with a 

particular solution, all three components of validation uncertainty would need to be evaluated.

5.1.2 Code Verification

Code verification evaluates whether the mathematical formulation and constitutive relationship 

are correctly coded in the program. Code developers generally perform this exercise. One way 

of verifying that coding is correct is by simulating a problem that has an exact solution and

evaluating whether the code reproduces the exact solution. Generally, code verification should 

precede solution verification to ensure that an error-free code is being used to obtain the 

solution of interest. In the present analysis, a commercial CFD package was used and 

therefore code verification was not performed.

5.1.3 Solution Verification

According to ASME V&V 20-2009 [ASME 2009], solution verification refers to the process of 

estimating the numerical uncertainty for a particular solution of a problem. The two main 

sources of errors that contribute to numerical uncertainty are the discretization error and the 

iteration convergence error. The discretization error is the difference between the result of a 

particular simulation with a specific grid and the results obtained with systematically refined 

grids. Discretization error can be estimated using a systematic grid refinement analysis, where 

the numerical solution is expected to asymptotically converge to the exact solution as the mesh 

size approaches zero. Systematic grid refinement forms the basis of the solution 

verification process.

Most CFD codes use iterative methods in which the numerical result must converge to the exact 

value as the iterations develop. Iteration convergence is related to the number of iterations 

required to obtain residuals that are sufficiently close to zero, either for a steady-state problem 

or for each time step in an unsteady problem. Iteration error is estimated using the residual root 

mean square (RMS) between subsequent iterations of a variable over all the volumes of the 

computational domain. In general, solutions are checked for iteration convergence by 

monitoring critical variables as solutions iteratively progress. It is necessary to achieve 

2

exp

2

input

2

uval = unum + u + u
Uncertainty Quantification and Grid Convergence Index (GCI) Analysis

5-3

complete iteration convergence before any discretization error estimate is performed because 

incomplete iteration convergence will adversely affect the estimate [Roache 1998, 1997, 1994]. 

Discretization error is normally calculated using Richardson extrapolation (RE) 

[Richardson 1910, Richardson and Gaunt 1927]. Richardson extrapolation is formulated to 

estimate the exact solution to the partial differential equations (PDEs) using the order of 

accuracy and solutions from two levels of grids. Richardson extrapolation is given by

φh = φexact + Chp Equation (5-2)

where 

φh = Numerical solution 

φexact = Exact solution to the PDEs

h = Grid size

C = Convergence coefficient 

p = Order of accuracy 

and the discretization error Edis is given by

Edis = (h) – exact = Chp Equation (5-3)

Using Equation (5-3), it is possible to determine the discretization error using solutions from two 

grid levels if the order of accuracy is known. For most practical problems, however, the order of 

accuracy realized from the solution is not the same as the formal order of accuracy specified in 

the model and is not known in advance. Hence, solutions from at least three levels of mesh are 

needed to determine the actual order of accuracy, p, and the discretization error. Moreover, 

these grids have to be sufficiently fine such that the grid is in the asymptotic range to obtain a 

reliable value of p. It is difficult to generate grids in the asymptotic range for complex 

engineering problems. Researchers have found a number of disadvantages in directly using RE 

to calculate numerical uncertainty [Roache 1998]. Based on the concept of RE, Roache 

proposed a method for a numerical uncertainty estimate, known as the grid convergence index 

(GCI) method [Roache 1994, 2003, 1998; Pelletier and Roache 2006]. The GCI method was 

developed as a means of converting the discretization error estimates from the RE into 

uncertainties. The GCI is an estimated 95 percent uncertainty calculated by multiplying the 

absolute value of the RE error by a factor of safety (Fs) [ASME 2009]. GCI provides the 

uncertainty of the solutions derived from discretized equations by a numerical and an analytical 

approach. The GCI calculation formula is detailed in subsequent sections. The present study

uses the GCI technique to estimate the numerical uncertainty for the temperature values 

calculated in HSM–1 and HSM–15.

5.2 Grid Convergence Index Calculation Method 

This section describes the sequential steps needed to calculate GCI for three levels of grids. 

Celik et al. [2008] and ASME V&V 20-2009 [ASME 2009] described the procedure in detail. The 

procedure involves the following steps.

Step 1: Define a representative cell, mesh, or grid size h. It can be calculated as

Equation (5-4)

φ φ

1/ 3 N

i 1

h Vi / N



 



 



 



 = ∑∆

=
Uncertainty Quantification and Grid Convergence Index (GCI) Analysis

5-4

where 

N = Total number of grids/cells used in the computations

ΔVi = Volume of the ith cell

Step 2: Generate three significantly different grids through systematic grid refinement. It is 

desirable that the mesh refinement factor (r), defined as hcoarse/hfine, be greater than 1.3. The 

grid refinement should be done in a systematic way (i.e., the refinement itself should be 

structured even if the grid is unstructured). Once the grids have been generated, simulations 

should be run to determine the values of key variables important to the objective of the 

GCI analysis.

Step 3: For the three grids with representative cell sizes h1, h2, and h3, with h1 < h2 < h3, define 

the mesh refinement factors r21 = h2/h1, r32 = h3/h2 and calculate the apparent (or observed) 

order, p, of the method using the following three equations

[ ( )] 



 



 +

ε

ε

p = 1/ln r ln q(p)

21

32

21

Equation (5-5)

Equation (5-6)

Equation (5-7)

where ε32 = φ3- φ2 ε21 = φ2- φ1 and φk denotes the simulation value of the variable on the kth

grid level.

The parameter “s” in Equation (5-6) is the indicator of oscillatory convergence in the solution. 

When s = 1, the grid convergence under consideration is monotonic. When s = −1, the grid 

convergence is oscillatory. By considering “s” in Equation (5-6), one ensures that the calculation 

procedure provides GCI even if the grid convergence under consideration is oscillatory 

in nature.

Equations (5-5) through (5-7) can be solved using fixed point iteration with the initial guess 

equal to the first term (i.e., q = 0) [ASME 2009].

Step 4: Calculate the extrapolated values from the following equations

Equation (5-8)

Equation (5-9)

Step 5: Calculate the following error estimates along with the observed order of the method p.

 Approximate relative error:

















−

− = r s

r s

q(p) ln

p

32

p

21





 





ε

ε = • 21

32 s 1 sign

( )

(r 1)

r

p

21

1 2

p

21 21

ext −

φ − φ φ =

( )

(r 1)

r

p

32

2 3

p

32 32

ext −

φ − φ φ =
Uncertainty Quantification and Grid Convergence Index (GCI) Analysis

5-5

Equation (5-10)

Equation (5-11)

Extrapolated relative error: 

Equation (5-12)

The fine grid convergence index is calculated as

Equation (5-13)

In the previous expression, Fs is the factor of safety that is introduced to ensure a GCI with a 

95 percent confidence interval. Originally, Fs was assigned a value of 3 for two-grid studies 

[Roache 1994], but Roache [1998] recommended a less conservative value of 1.25 for cases 

using at least three grid solutions and the observed p. According to ASME V&V 20-2009

[ASME 2009], a value of 1.25 should be used for Fs with three-grid studies involving structured 

grid refinement. 

Celik et al. [2008] modified Equation (5-5) in their GCI formula, given by

[ ( )] 



 





+

ε

ε

p = 1/ln r ln q(p)

21

32

21

Equation (5-14)

According to Celik and Karatekin [1997], the absolute value of ε32/ε31 is necessary to ensure 

extrapolation toward h = 0.0. Negative values of ε32/ε31 < 0.0 indicate oscillatory convergence.

5.3 Uncertainty Quantification and GCI Analysis for HSM–1 Module

Solution verification for both HSM–1 and HSM–15 modules was performed using four levels of 

mesh generated with systematic grid refinement. Simulations were performed using all four 

levels of mesh at three specific time points, which are calendar years 1993, 2012, and 2093. 

The objective of this analysis was to determine the numerical uncertainty due to discretization 

error based on the GCI technique with a 95 percent confidence interval containing the exact 

solution of the problem and system of equations. Simulations were performed for three 

calendar years to assess the variation of the apparent order p and GCI values with decreasing 

thermal load. The levels of grids used in the analysis for the HSM–1 module are provided in 

Table 5-1, which also lists the mesh refinement factor between two consecutive levels of grids

calculated using the procedure described in Step 1 and Step 2 in Section 5-2. For example, the 

mesh refinement factor r for the baseline mesh is calculated using the grid sizes of the baseline 

and coarse meshes. Simulation for all the levels of grids was performed with the same input 

2

32 2 3

a e

φ

φ − φ =

1

21 1 2

a e

φ

φ − φ =

12

ext

1

12

21 ext

ext e

φ

φ − φ =

r 1

F e GCI

p

21

21

21 s a

fine −

• =
Uncertainty Quantification and Grid Convergence Index (GCI) Analysis

5-6

parameters and the same turbulence model that was used for the baseline simulations for the 

HSM–1 module described in Chapter 2 of this report.

Table 5-1. Mesh levels used in GCI analysis of HSM–1 module

Mesh type Total cell number Mesh refinement factor (r)

Coarse mesh 1,220,596

Baseline mesh 3,528,864 1.4245

Fine mesh 1 9,214,367 1.3771

Fine mesh 2 25,000,000 1.3856

Figure 5-1 shows temperature contours along the transverse mid-plane for the four mesh levels 

for HSM–1 simulations at the representative year 1993. It can be observed that there is no 

significant qualitative difference in the simulated results from the four mesh levels, indicating 

that mesh refinement does not cause any major deviation in the simulated temperature field. 

However, temperature contours show that the finest mesh calculation captures the thermal 

gradient better compared to the coarse mesh and baseline mesh simulations. 

Tables 5-2, 5-3, and 5-4 provide the simulated values of the maximum and minimum cladding 

temperatures for the HSM–1 module for years 1993, 2012, and 2093, respectively, for all four 

mesh levels. There are only minor differences in the predicted values of the maximum and 

minimum cladding temperatures from the different grid levels for the HSM–1 module. The 

maximum temperature difference between the coarse mesh solution and the finest mesh 

solution is approximately 3 K [5.4 °F] and in most of the cases the variation is less than 

1 K [1.8 °F]. Also note that the variation decreases with increased storage time.

To better understand the variability of predicted temperature values, the percentage of variation 

of the maximum and minimum cladding temperatures with mesh size is plotted in Figure 5-2 for 

the different grid levels for the HSM–1 module. The percentage of variation is calculated using 

the following relationship

Variation (∆V) = ( ) x 100

Equation (5-15)

The x-axis in Figure 5-2 shows the nondimensionalized mesh size (h/hfine-mesh-2), where h is the 

representative grid size as defined in Step 1 and hfine-mesh-2 is the corresponding grid size for the

fine mesh 2 case. Hence, a nondimensional mesh size of 1 corresponds to the finest mesh 

(fine mesh 2), with higher values for coarser meshes. The calculated deviation is shown along 

the y-axis. Note that the variation between the solutions from different grid levels is within 

1 percent.

The apparent order and the GCI values are calculated for the variables maximum and minimum 

cladding temperatures. Equation (5-13) was used to calculate the GCI using the modified 

apparent order p described in Equation (5-14). For the GCI calculations, the factor of safety Fs

is taken as 3.0. The values for p and the GCI values for the HSM–1 modules are provided in 

Table 5-5. The GCI values listed in the table were calculated based on the baseline and two 

fine-mesh-level solutions. 

finest mesh

finest mesh

T

T T

−

− −
Uncertainty Quantification and Grid Convergence Index (GCI) Analysis

5-7

(a) Coarse mesh temperature 

contours

(b) Baseline mesh temperature contours

(c) Fine mesh 1 temperature contours (d) Fine mesh 2 temperature contours

Figure 5-1. Temperature contours along the mid-plane for four mesh levels for 

simulations of HSM–1 at year 1993 
Uncertainty Quantification and Grid Convergence Index (GCI) Analysis

5-8

Table 5-2. Computed maximum and minimum cladding temperature for 

HSM–1 module (1993)

Mesh type

Maximum cladding temperature 

(K) °F = [1.8 K-460]

Minimum cladding temperature (K)

°F = [1.8 K-460]

Coarse mesh 421.661 319.41

Baseline mesh 420.363 318.299

Fine mesh 1 421.518 317.69

Fine mesh 2 419.399 316.38

Table 5-3. Computed maximum and minimum cladding temperature for 

HSM–1 module (2012)

Mesh type

Maximum cladding temperature

(K) °F = [1.8 K-460]

Minimum cladding temperature

(K) °F = [1.8 K-460]

Coarse mesh 392.017 315.595

Baseline mesh 390.82 314.07

Fine mesh 1 391.676 314.442

Fine mesh 2 390.963 313.577

Figure 5-2. Deviation of temperature values calculated at different mesh levels for 

HSM–1: deviations are calculated as percentage of fine mesh temperatures

Table 5-4. Computed maximum and minimum cladding temperature for 

HSM-1 module (2093)

Mesh type

Maximum cladding temperature

(K) °F = [1.8 K-460]

Minimum cladding temperature

(K) °F = [1.8 K-460]

Coarse mesh 343.209 308.437

Baseline mesh 342.786 307.659

Fine mesh 1 342.954 307.539

Fine mesh 2 342.537 307.254
Uncertainty Quantification and Grid Convergence Index (GCI) Analysis

5-9

Table 5-5. Apparent order (p) and calculated GCI for HSM–1 module

Year Variable Apparent order (p) GCI

1993

Maximum cladding temperature 1.83 1.85

Minimum cladding Temperature 2.26 1.14

2012

Maximum cladding temperature 0.56 2.67

Minimum cladding Temperature 2.55 0.84

2093

Maximum cladding temperature 2.75 0.25

Minimum cladding Temperature 2.56 0.21

The values of p and GCI provide an evaluation of numerical uncertainty in the simulations. The 

values of GCI (less than 4 percent in most of the cases) in Table 5-5 indicate that there was no 

significant deviation of the calculated results due to mesh refinement or mesh convergence, 

further strengthening the observations made in Figures 5-1 and 5-2. For all the simulations, the 

spatial discretization used was second-order accurate. It can be observed from Table 5-5 that 

the apparent order of accuracy is not constant and is, in most of the cases, higher than the 

theoretical order of accuracy (p = 2).

Figure 5-3 shows the variation of calculated p values for the maximum and minimum cladding 

temperatures at the three calendar years plotted against the nondimensional mesh size 

(h/hfine-mesh-2). The apparent order at a mesh point was calculated using a combination of three 

temperatures. The following combinations were used for the apparent order calculations: 

(i) Coarse, baseline, and fine mesh 1 solutions (mesh set 1).

(ii) Coarse, fine mesh 1, and fine mesh 2 solutions (mesh set 2).

(iii) Baseline, fine mesh 1, and fine mesh 2 solutions (mesh set 3).

It is possible to have other combinations of solutions, but these three sets provide 

representative values as they generally consider progressively refined meshes. The purpose of 

this plot (Figure 5-3) was to understand the general variability of p with mesh refinement. From 

the figure, it is clear that the calculated value of p varies significantly for different levels of mesh 

refinement.

From Figure 5-3 it can also be observed that the calculated values of p are significantly higher in

many cases. The theoretical value of p is shown by the red dotted line in Figure 5-3. The 

dispersivity of the calculated p values is significant for the current simulations. According to 

Roache [1998], the general disagreement of the calculated p-value and the theoretical p-value is 

not necessarily a sign that the simulated values or simulations are wrong and unsatisfactory; 

however, further analysis was performed to understand this behavior. 

5.3.1 Assessment of Asymptotic Behavior of the Solution

To find a probable explanation for this behavior, two quantities were calculated and their 

variations with mesh size were analyzed. The first quantity under consideration is the 

discretization error that is defined in Equation (5-3).
Uncertainty Quantification and Grid Convergence Index (GCI) Analysis

5-10

Figure 5-3. Calculated apparent order using different mesh levels for HSM–1

Using Equation (5-3), the convergence coefficient can be defined as 

C = Edis/hp Equation (5-16)

As described in ASME V&V 20-2009 [ASME 2009], for a well-behaved problem the values of 

the convergence coefficient c should become constant as the grid is refined for a uniformly pth

order method. Hence, the slope of the Edis distribution against h should be continuous 

and constant. 

An exact analytical solution for the dry cask thermal problem, however, is not known; according 

to Roache [2003], the extrapolated solution can be used to calculate the discretization error 

when the exact solution is not known. In the present study, the discretization error is computed 

using the extrapolated values obtained from two different solutions using Equation (5-17).

Equation (5-17)

For the present study, therefore, the discretization error is calculated using Equations (5-3) and 

(5-17). From these equations, it is clear that to calculate the discretization error, three sets of 

solutions are needed. 

Figure 5-4 shows the variation of the discretization error with nondimensional mesh size 

(1 corresponds to finest mesh) at three calendar years (1993, 2012, and 2093). The 

discretization error at every point was calculated using a set of three solutions. The same 

combination of solutions described in connection with Figure 5-3 was used to determine 

discretization errors for Figure 5-4. Note that the curves do not exhibit asymptotic behavior and 

have sharp discontinuities. To further illustrate the nonasymptotic nature of the discretization 

r 1

r

p

21

1 2

p

21

ext

−

φ − φ φ =
Uncertainty Quantification and Grid Convergence Index (GCI) Analysis

5-11

Figure 5-4. Distribution of discretization error with mesh refinement for HSM–1

error and the solution, the variation of the convergence coefficient [C in Equation (5-16)] due to 

mesh refinement is shown in Figure 5-5. The convergence coefficients were calculated using 

the same combination of mesh and solutions that were used to determine discretization errors 

for Figure 5-4. If the solutions are within the asymptotic region, then the convergence

coefficient C [slope of Equation (5-3)] should be a constant for all mesh sizes. However, one 

can see that the variation of C doesn’t follow any general pattern. The value of C varies 

significantly, which indicates oscillatory convergence. Figures 5-4 and 5-5 suggest that the 

solutions obtained from the systematic mesh convergence study may not all lie in the 

asymptotic convergence region.

5.3.2 Iteration Convergence

As indicated in Section 5.1.3, one of the important factors in numerical uncertainty estimation 

calculation is iteration convergence. The possibility of incomplete iteration convergence was 

explored for the present study. Figure 5-6 shows the variation of the average temperature of

selected monitor points with iteration for the baseline case at calendar year 2012. The plot

shows an overall convergence with iteration, and the trend indicates that the solution is 

converged. Figure 5-7 highlights the convergence pattern for the last 6,000 iterations as 

function of temperature and indicates some minor oscillations. Although one can conclude from 

Figure 5-6 that the solution has an overall convergence, Figure 5-7 shows that some minor 

oscillations may still persist. Hence, the calculation of p and GCI essentially depends on the

iteration at which the solution was terminated. 

5.3.3 Other Numerical Issues

It is clear from Tables 5-2, 5-3, and 5-4 that the calculated temperature values do not change 

significantly due to mesh refinement and, generally the variation is less than 1K [1.8 0

F].

Essentially, these differences can be attributed to numerical noise and machine error. Because 
Uncertainty Quantification and Grid Convergence Index (GCI) Analysis

5-12

Figure 5-5. Variation of convergence coefficient with mesh refinement for HSM–1

Figure 5-6. Variation of average temperature of selected monitor points with iteration,

showing overall convergence for HSM-1
Uncertainty Quantification and Grid Convergence Index (GCI) Analysis

5-13

Figure 5-7. Variation of average temperature of selected monitor points with iteration for 

last 6,000 iterations, highlighting minor oscillations for HSM-1

the values predicted by different mesh levels are very close to each other, such minor 

oscillations tend to affect some of the terms in Equations 5-8 through 5-14 disproportionally and 

subsequently affect the calculation of apparent order p and GCI values.

According to Roache [1998], it is difficult to have complete iteration convergence in multiphysics 

problems with multiple scales of solution and multiple scales in turbulence modeling because 

different parts of the physics may converge differently. Roache [1998, 1997] indicates that, in 

multiphysics problems, a grid may be in the asymptotic range for some flow physics, but not in 

the asymptotic range for other physics. This behavior can be more prominent in problems that 

involve the laminar to turbulent transition. The present problem involves flow that may be in the 

transition region, which contributes to the uncertainty. In addition, the present problem involves 

a complicated geometry that influences the flow path and causes repeated flow separations and 

reattachments that are difficult to model accurately [Roache 1997]. Additionally, an asymptotic 

solution requires a systematic mesh refinement in all three coordinate directions. In the present 

study, an attempt has been made to achieve a systematic mesh refinement in all directions with 

the same mesh refinement factor. However, due to the complexity of the geometry and certain 

limitations of the grid generation package, the mesh refinement factor had to be adjusted in 

small parts of the domain—especially in narrow gaps and annuli. Therefore, the complexity of 

the flow physics and limitations in mesh refinement may have resulted in nonasymptotic 

solutions. However, the closeness of the calculated temperature values showed that mesh 

refinement may have had little impact on the solution.
Uncertainty Quantification and Grid Convergence Index (GCI) Analysis

5-14

5.3.4 Modified Estimation of GCI for HSM-1

To obtain a useful bound on the simulated results and a realistic numerical uncertainty estimate, 

a modified approach to calculate GCI was used. In this approach, a theoretical order of 

accuracy p=1 was used instead of the apparent order of accuracy. This is because the 

simulations involved specification of mixed order methods, whereby the spatial discretization 

was performed using a second order scheme and the body force averaging in the pressure 

equations was performed using a first order method. In addition, the radiative transport 

equations in the radiation discrete ordinate model were solved using a first order method. 

Therefore, the overall solution would be contaminated by these first order approximations and 

not be second order accurate. Hence, the theoretical order of accuracy for the GCI calculation 

was assumed to be 1.

The general objective of the GCI calculations described in this chapter was to provide an idea of 

the numerical uncertainty associated with the benchmark modeling (Chapter 3) and the thermal 

evolution study over 300 years (Chapter 4). As these studies were performed using the 

baseline mesh, it was imperative to use the baseline mesh results for the GCI calculations. 

Therefore, Equations (5-11) and (5-13) were used to calculate the GCI with p=1. The mesh 

refinement factor r21 and the relative error term ea

21 were calculated based on baseline and 

coarse mesh simulations. The calculated values of GCI for HSM-1 using this approach are 

provided in Table 5-6.

The calculated GCI values are greater for this bounded case compared to the values provided 

in Table 5-5 but are still low, indicating the minor influence of mesh refinement.

5.4 Uncertainty Quantification and GCI Analysis for HSM–15 Module 

The methodology used to perform uncertainty quantification and GCI calculation for HSM-1 was 

applied also to HSM-15. The four levels of meshes used in the analysis for HSM–15 and the 

mesh refinement factors are provided in Table 5-7. The mesh refinement factors were 

calculated using steps 1 and 2 described in Section 5-2 using two consecutive levels of meshes. 

All four levels of meshes were simulated with the same input parameters and the same 

turbulence model that were used for the baseline simulations for the HSM–15 module. As for

HSM-1, simulations for HSM-15 at every mesh level were performed at three different calendar 

years, which are years 1996, 2012, and 2097. 

Tables 5-8, 5-9, and 5-10 provide the simulated values of the maximum and minimum cladding 

temperatures for the HSM–15 module for years 1996, 2012, and 2097, respectively. One can 

generally observe that, similar to the case for the HSM–1 module, there are only minor 

differences in the predicted values of the maximum and minimum cladding temperatures for the 

different grid levels for the HSM–15 module. The largest temperature difference between the 

coarse and fine mesh solutions, observed for the minimum cladding temperature at year 2012, 

is less than 3 K [5.4 °F]. In most of the cases the difference is less than 1 K [1.8 °F]. As pointed 

out in Sections 5.3.2 and 5.3.3, these minor differences could be artifacts of machine round-off 

errors and oscillatory iteration convergence. Figure 5-8 shows the percentage variation of the 

maximum and minimum cladding temperatures calculated using Equation (5-15) between 

solutions for the different grid levels for HSM–15. As can be observed, there is insignificant 

variation of less than 1 percent between the solutions with the different grid levels. While the 

variability decreases with mesh refinement, the close similarity of computed results shows that 

the effect of mesh refinement was not significant.
Uncertainty Quantification and Grid Convergence Index (GCI) Analysis

5-15

Table 5-6. GCI value for a conservative bounded value of p = 1 for HSM–15 module

Year Variable

Value from baseline 

mesh simulation 

(K) °F = [1.8 K-460] GCI

1993

Maximum Cladding temperature 420.363 2.181

Minimum Cladding Temperature 318.299 2.466

2012

Maximum Cladding temperature 390.82 2.16

Minimum Cladding Temperature 314.07 3.43

2093

Maximum Cladding temperature 342.786 0.872

Minimum Cladding Temperature 307.659 0.178

Table 5-7. Mesh levels used in GCI analysis of HSM–15 module

Mesh Type Total cell number Mesh refinement factor (r)

Coarse mesh 1,136,626

Baseline mesh 3,278,964 1.4235

Fine mesh 1 9,000,000 1.3744

Fine mesh 2 22,635,832 1.3853

Table 5-8. Computed maximum and minimum cladding temperature for 

HSM–15 module (1996)

Mesh type

Maximum cladding temperature

(K) °F = [1.8 K-460]

Minimum cladding temperature 

(K) °F = [1.8 K-460]

Coarse mesh 502.23 327.56

Baseline mesh 502.61 326.49

Fine mesh 1 503.04 325.19

Fine mesh 2 502.4 324.96

Table 5-9. Computed maximum and minimum cladding temperature for 

HSM–15 module (2012)

Mesh type

Maximum cladding temperature

(K) °F = [1.8 K-460]

Minimum cladding temperature

(K) °F = [1.8 K-460]

Coarse mesh 463.13 323.29

Baseline mesh 462.83 321.66

Fine mesh 1 462.21 321.40

Fine mesh 2 461.63 320.50

Table 5-10. Computed maximum and minimum cladding temperature for 

HSM–15 module (2097)

Mesh type

Maximum cladding temperature

(K) °F = [1.8 K-460]

Minimum cladding temperature

(K) °F = [1.8 K-460]

Coarse mesh 364.77 310.22

Baseline mesh 364.70 309.21

Fine mesh 1 364.20 309.56

Fine mesh 2 364.08 308.87
Uncertainty Quantification and Grid Convergence Index (GCI) Analysis

5-16

Figure 5-8. Deviation of temperature values calculated at different mesh levels for 

HSM–15: deviations are calculated as percentage of fine mesh temperatures

The values for apparent order p and GCI for HSM–15, calculated using the baseline, fine mesh 

1, and fine mesh 2 solutions, are provided in Table 5-11. The values of calculated GCI indicate 

again that the variability between the solutions obtained at different mesh levels. As observed 

in connection with HSM-1, the apparent order of accuracy p for HSM-15 is not constant and 

some of the calculated values are higher than the theoretical order of accuracy. Figure 5-9 

shows the calculated p values for the maximum and minimum cladding temperatures at the 

three calendar years for HSM-15. As was the case for HSM-1, the apparent order at a mesh 

point was calculated using a combination of three temperatures. The following combinations 

were used for the apparent order calculations: 

(i) Coarse, baseline, and fine mesh 1 solutions (mesh set 1).

(ii) Coarse, fine mesh 1, and fine mesh 2 solutions(mesh set 2).

(iii) Baseline, fine mesh 1, and fine mesh 2 solutions (mesh set 3).

The general objective of this plot was to determine whether the value of p remains constant and 

lies in the asymptotic region. The figure indicates that p is not constant and varies significantly.

Further analysis was performed to understand solution behavior variation with mesh refinement. 

Figure 5-10 shows the distribution of the discretization error for different mesh sizes for HSM–15

at three calendar years. The discretization error was computed using Equation (5-3) and the 
Uncertainty Quantification and Grid Convergence Index (GCI) Analysis

5-17

Table 5-11. Apparent order (p) and calculated GCI for HSM–15 module

Year Variable Apparent order GCI

1996 Maximum cladding temperature 1.2 0.79

Minimum cladding temperature 5.48 0.04

2012 Maximum cladding temperature 1.94 0.423

Minimum cladding temperature 3.6 0.37

2097 Maximum cladding temperature 4.5 0.03

Minimum cladding temperature 2.05 0.70

Figure 5-9. Calculated apparent order using different mesh levels for HSM–15

Figure 5-10. Distribution of discretization error with mesh refinement for HSM–15
Uncertainty Quantification and Grid Convergence Index (GCI) Analysis

5-18

same combination of meshes used in Figure 5-9. It is clear from the figure that the behavior is 

not asymptotic with mesh refinement. The convergence coefficient C is plotted in Figure 5-11 

using the same combination of meshes. The figure shows that the pattern of distribution for C is 

also not asymptotic. Hence, Figures 5-10 and 5-11 are suggestive of oscillatory convergence. 

The possible causes of this pattern are similar to those observed in connection with HSM–1, as 

discussed in Sections 5.2.2 and 5.2.3. These potential explanations include incomplete iteration 

convergence, numerical effects such as machine error, complex geometry and flow physics, 

and limitations in systematic mesh refinement. Since the internal geometrical configurations of 

HSM-1 and HSM-15 are the same and both involve the same basic flow physics, these factors 

had similar effects on the HSM-15 solutions and contributed to the nonasymptotic patterns with 

mesh refinement.

As for HSM-1, a modified approach to calculate GCI for HSM-15 was adopted in order to obtain 

a useful bound of the simulated results and a realistic numerical uncertainty estimate. The 

detailed rationale behind the modified approach was explained in Section 5.3.4. A theoretical 

order of accuracy p=1 was used, instead of the apparent order of accuracy, in order to reflect 

the fact that some of the equations were solved using first order approximations that would 

affect the overall accuracy of the solution. In addition, the baseline mesh was used to calculate 

the GCI values because the benchmark modeling, parametric analysis, and thermal evolution 

study used the baseline mesh. The mesh refinement factor r21 and the relative error term ea

21

were calculated using the baseline and coarse mesh simulations. The calculated values of GCI 

for HSM-15 using this approach are provided in Table 5-12. Again, the small calculated GCI 

values, even using a lower order of accuracy, indicate limited variability due to mesh refinement.

5.5 Observations from GCI Analysis

The major observations from the GCI and uncertainty analysis are as follows.

• For maximum and minimum cladding temperatures, differences among the solutions 

from different mesh levels are not significant.

• The apparent order of accuracy calculated using different combinations of solutions was 

not constant and was in many of the cases higher than the theoretical order of accuracy

(p = 2).

• A detailed discretization convergence analysis showed that the solutions may not all lie 

in the asymptotic convergence zone.

• A number of factors may have influenced the calculated order of convergence p and 

GCI. These include (i) minor variation in the solution even when the system achieved an 

overall iteration convergence and (ii) machine errors.

• The present problem involves turbulent natural convection through a geometrically 

intricate domain, in which different aspects of the physics may converge differently,

making the determination of GCI a challenging exercise. As mentioned in Section 2.5, 

the flow behaviors in HSM–1 and HSM–15 both involve the laminar-turbulent transition. 

With the higher heat load in HSM–15, the flow is even less likely to be in the fully 

turbulent regime. A modified approach, with a low apparent order of 1 and a high factor 

of safety, was employed to calculate the GCI.
Uncertainty Quantification and Grid Convergence Index (GCI) Analysis

5-19

• The calculated GCI values showed minor variation in results, indicating that mesh 

refinement may not have a significant impact on temperature calculations for HSM–1 

and HSM–15.

Figure 5-11. Variation of convergence coefficient with mesh refinement for HSM–15

Table 5-12. GCI value for a conservative bounded value of p = 1 for HSM–15 module

Year Variable

Value from baseline 

mesh simulation 

(K) °F = [1.8 K-460] GCI

1996

Maximum cladding temperature 502.61 0.534

Minimum cladding temperature 326.49 2.314

2012

Maximum cladding temperature 462.83 0.13

Minimum cladding temperature 321.66 2.33

2097

Maximum cladding temperature 364.70 0.469

Minimum cladding temperature 309.21 3.59

6-1

6. SUMMARY AND CONCLUSIONS

A three-dimensional (3-D) computational fluid dynamics (CFD) model was developed for 

ventilated horizontal storage modules HSM–1 and HSM–15 at the Calvert Cliffs Independent 

Spent Fuel Storage Installation (ISFSI). These horizontal storage modules contain a modified 

version of the NUHOMS dry storage canister (DSC) that stores 24 pressurized water reactor 

(PWR) fuel assemblies. Each assembly accommodates a 14 × 14 arrangement of Combustion 

Engineering (CE) fuel rods. The CFD model, developed using the commercial package 

ANSYS-FLUENT® Version 14.5, explicitly models the concrete storage units, DSC shells, and 

internal structures such as spacer plates, tie rods, fuel assemblies, and the supporting rails. 

The storage system was meshed using unstructured hexahedral grids with mesh refinement in 

regions expected to have higher gradients of temperature and velocity. Individual fuel rods 

within a fuel assembly were not represented explicitly; the volume within each fuel assembly

was modeled as a continuous porous medium. A separate 3-D CFD model of a fuel assembly

was developed to determine the inertial and frictional flow resistance factors, which were used 

in the study to account for the pressure drop within the fuel assemblies. A separate calculation 

based on the CE 14 × 14 fuel configurations was done to determine the orthotropictemperature-dependent thermal conductivity in the fuel assembly area. For the baseline case 

simulations, the low Reynolds number k-ε model with buoyancy effect was used for the 

turbulence model. The discrete ordinate (DO) radiation model was used for radiation. The 

decay heat was specified as a volumetric heat generation rate within the fuel assemblies and 

was applied over the active fuel length within each fuel region and modified using an axial 

peaking factor to include axial variation.

6.1 Baseline Simulations and Comparison of Simulated to 

Measured Temperatures

Baseline simulations for HSM–1 and HSM–15 were performed to understand the flow and 

temperature field within the storage system. Simulation results yielded buoyancy-driven natural 

convection flow patterns within the module, in which cooler air enters the system through the 

opening and absorbs the decay heat, lowering the air density, and the warmed air exits through 

the outlet vent at the top of the structure. The simulated patterns of temperature distribution for 

the DSC shell are similar for HSM–1 and HSM–15. The temperature gradually increases 

toward the mid-section of the DSC from relatively low values at the front and back end sections. 

This pattern can be attributed to the decay heat load distribution and axial peaking factor and to 

the asymmetric airlow around the DSC in the axial direction. The high-temperature region on 

the upper half of the curved surface of the DSC shell is due to the buoyancy-driven upward flow 

of heated air around the DSC and helium flow within the DSC. The highest simulated 

temperature for the HSM–15 module DSC shell {391.46 K [245.2 °F]} is significantly higher than 

the highest simulated temperature for the HSM–1 DSC shell {355 K [180 °F]}. This is principally 

because the decay heat load is higher in the HSM–15 module. Another factor likely contributing 

to the lower overall temperature in HSM–1 is the fact that one side of the HSM–1 module is 

exposed to the atmosphere.

Computed temperatures for the HSM–1 and HSM–15 modules were compared to temperatures 

measured on June 27 and 28, 2012, as part of the Used Fuel Disposition Campaign of the 

U.S. Department of Energy [Suffield et al. 2012]. Due to constraints on access to the different 

regions of the DSC and safety constraints, temperature measurements were obtained only on 

the DSC shell outer surfaces. The front door of each module was opened to obtain the limited 

temperature measurements. The simulated temperatures for both modules were, on average, 

10–12 K [18–22 °F] higher than the measured values. The differences between the computed 
Summary and Conclusions

6-2

and measured values were smaller on the end cap of the DSC shell and the differences were

larger for locations away from the end cap. The measured temperatures have significant 

uncertainty because of the difficulties in obtaining accurate measurements in the narrow gap 

between the DSC and storage module wall and may reflect mixing of ambient air with the 

internal air through the open doors. The measurements also show a relatively flat temperature 

distribution along the axial direction of the DSC shell top wall; this distribution is unlikely given 

the magnitude and axial distribution pattern of the heat load. Hence, a thorough assessment of 

experimental uncertainty is needed to make a meaningful comparison with computed data and 

use it to assess the validity and accuracy of the numerical solution.

6.2 Sensitivity Analyses

A series of sensitivity studies was performed using the baseline simulation. The following 

features were studied: (i) near wall mesh refinement, (ii) turbulence model, (iii) porous media 

resistance, (iv) insolation loading, and (v) extension of computational domain.

A study was conducted by adapting and refining the mesh near the DSC shell outer wall to 

understand the effects of y+ and mesh spacing on temperature predictions. Three levels of 

mesh refinement were performed, with surface averaged y+ of 2.974, 1.575, and 0.8193. A 

comparison of temperature distributions obtained from simulations using these different meshes 

showed nearly identical values and patterns. Minor deviations were observed near the top of 

the DSC shell, where the flow is characterized by higher gradients. In general, the close 

agreement among the simulated temperatures demonstrates that the baseline mesh is 

adequate for the present study.

The turbulence models considered in the sensitivity analysis included the baseline low 

Reynolds number k-ε model, the standard k-ε model with scalable wall function and buoyancy 

effects, the shear stress transport k-ω model with low Reynolds number corrections, and the 

shear stress transport k-ω model with buoyancy effects. Sensitivity analyses for the turbulence 

models showed that the simulated results with the baseline low Reynolds k-ε model predicted 

the lowest temperatures. 

For the sensitivity analysis on the porous media model, four cases were simulated: higher 

porous media resistance than the baseline, lower porous media resistance than the baseline, 

zero resistance, and no porous media assumption. The study revealed that varying porous 

media resistance to a higher or lower value compared to baseline did not affect the temperature 

distribution. However, studies with zero resistance in the porous media and with no porous 

media showed a relatively flat temperature distribution pattern that is markedly different from the 

other cases. 

A sensitivity analysis was performed to understand the effect of insolation. Comparing

temperatures from the baseline simulations (which consider insolation) with temperatures from 

simulations that omit insolation, the HSM component temperatures specified on the outer 

surfaces are sensitive to insolation by several degrees, although extreme temperatures are 

less sensitive.

A separate study was performed to check the accuracy of the convective heat transfer 

coefficient specified as the boundary condition on the outer walls of the storage modules. This 

was done by solving the baseline problem using an extended domain that included the 

surrounding atmosphere, where the external surface heat transfer was calculated from the flow 

variables in the solution without the need for an explicit specification of the boundary condition. 
Summary and Conclusions

6-3

The average convective heat transfer coefficient calculated from the extended domain study 

was comparable to the value used in the baseline simulations.

6.3 Simulation of Long-Term Storage Conditions

The baseline computational models for HSM–1 and HSM–15 were used to predict thermal 

behavior for 300 years of storage. Analyses for the 300-year duration used simulations 

assumed to be quasi-steady-state at multiple discrete points in time. The steady-state 

simulations were performed at 5-year intervals for the first 100 years and at 10-year intervals 

thereafter. 300-year simulations for both HSM–1 and HSM–15 configurations show (i) the drop 

in maximum cladding temperature with time is relatively rapid in the first 100 years and is

gradual thereafter and (ii) the minimum cladding temperature varies only a small amount as time 

progresses. Both these temperature distributions exhibited an asymptotic decay with increased 

time to a temperature slightly above the ambient temperature. The time evolution of 

temperature was also studied for other components, such as spacer plates, the fuel sleeve, and 

heat shield. All component temperatures showed similar trends of asymptotic decay to a 

temperature slightly above ambient, similar to the trend observed for cladding temperatures.

The general pattern of the temperature variations for the HSM–1 and the HSM–15 modules

were similar. However, the simulated temperatures for HSM–15 were consistently higher than 

the temperatures obtained from the HSM–1 simulations due to the higher decay heat load and 

location of the HSM–15 module.

6.4 Uncertainty Analysis

Finally, the numerical uncertainty of the simulated temperatures was evaluated using the grid 

convergence index (GCI) method. The GCI method is based on the principles of the classical 

Richardson extrapolation approach for calculating discretization error and is widely used to 

quantify the numerical uncertainty of the computed solution. In the present study the technique 

proposed by Roache [1998] and modified by Celik et al. [2008] was performed to determine the 

apparent order p and the GCI. The analysis used the extreme minimum and maximum 

temperatures for the cladding as convergence criteria of particular relevance to the performance 

of the storage module.

Four levels of mesh were used in the analysis: the baseline grid, a coarser grid obtained from 

systematically coarsening the baseline grid, and two levels of finer grids obtained by 

systematically refining the baseline grid. Simulations at three discrete time points within the 

300-year period were performed using these four grid levels for both HSM–1 and HSM–15 to 

obtain an accurate sense of the numerical uncertainties. The maximum cladding temperature 

and minimum cladding temperature were selected as the two representative variables for 

performing the GCI calculations. Other component temperatures are expected to have similar 

levels of numerical uncertainty as these two variables.

Computed temperature data obtained from simulation runs using different grid levels showed 

little variation. The apparent order of accuracy (p) calculated using different combinations of 

solutions at different mesh levels was not constant, indicating oscillatory convergence. In most 

of the cases, the calculated order p was higher than the theoretical order of accuracy. As a 

result, the GCI determined using these p values was low, indicating a small band of numerical 

uncertainty. This is most likely caused by the small differences in temperature from different 

mesh levels, which disproportionally affects some terms in the equation for apparent order of 

accuracy and results in an artificially high value of p. Other factors, such as minor oscillations in

a solution that shows overall convergence, may have affected the GCI calculations.
Summary and Conclusions

6-4

An alternative estimate of the numerical uncertainty was calculated using a theoretical order of 

accuracy p = 1 and a higher factor of safety of 3. This approach accommodated the fact that 

mixed order methods were used in some of the calculations and the radiative transport 

equations were solved using a first order technique. Because the benchmark study, parametric 

analysis, and thermal evolution study employed the baseline mesh, the final GCI calculation 

used the baseline mesh results. It was found that the GCI is still less than 4 percent, indicating 

that the level of discretization error in the solution is not significant and the associated numerical 

uncertainty of the results is low.
7-1

7. REFERENCES

Adkins, H.E., J.A. Fort, S.R. Suffield, J.M. Cuta, and B.A. Collins, “Thermal Modeling Studies 

for Active Storage Modules in the Calvert Cliffs ISFSI,” Proceedings of the International 

High-Level Radioactive Waste Management Conference April 28–May 2, 2013,

American Nuclear Society, Albuquerque, New Mexico, 2013.

ANSYS, Inc, “ANSYS FLUENT Version 14.5, User’s Guide,” ANSYS, Inc, 

Canonsburg, Pennsylvania, 2012.

ANSYS, Inc, “ANSYS CFDPOST Version 14.5, User’s Guide,” ANSYS, Inc. 

Canonsburg, Pennsylvania, 2012b.

Arya, M.S. and M. Keyhani, “Convective Heat Transfer in a Sealed Storage Cask Containing 

Spent-Fuel Canisters,” Nuclear Science and Engineering, Vol. 105, pp. 391–403, 1990.

ASME, “Standard for Verification and Validation in Computational Fluid Dynamics and Heat 

Transfer,” ASME V&V 20-2009, American Society of Mechanical Engineers, New York City, 

New York, 2009.

Bahney, R.H. and T.L. Lotz, “Spent Nuclear Fuel Effective Thermal Conductivity Report,” 

BBA000000-01717-5705-00010 REV00, Prepared by TRW Environmental Safety Systems, Inc., 

for the U.S. Department of Energy, Washington, DC, Jul 11, 1996.

Banken, G., K. Tavassoli, and J. Bondre, “Validation of Computational Fluid Dynamics Code 

Models for Used Fuel Dry Storage Systems,” Proceedings of the CFD for Nuclear Reactor 

Safety Applications (CFD4NRS-3) Workshop, Validation and Application of CFD and CMFD 

Codes to Nuclear Reactor Safety Issues, September 14–16, 2010, U.S. Nuclear Regulatory 

Commission, Washington, DC, 2010.

Bergman, T.L., A.S. Lavine, F.P. Incropera, and D.P. DeWitt, Fundamentals of Heat and Mass 

Transfer, Seventh Edition, John Wiley & Sons New York City, New York, 2011.

Canaan, R.E. and D.E. Klein, “A Numerical Investigation of Natural Convection Heat Transfer 

Within Horizontal Spent-Fuel Assemblies,” Nuclear Technology, Vol. 123, No. 2, pp. 193–208,

1998.

Celik, I. and O. Karatekin, “Numerical Experiments on Application of Richardson Extrapolation 

With Nonuniform Grids,” Journal of Fluids Engineers, Vol. 119, No 3, pp. 584–590, 1997.

Celik, I.B., U. Ghia, P.J. Roache, C.J. Freitas, H. Coleman, and P.E. Raad, “Procedure for 

Estimation and Reporting of Uncertainty Due to Discretization in CFD Applications,” ASME 

Journal of Fluids Engineering (Special Publication), Vol. 130, 2008.

Cha, J., S.N. Kim, and K.W. Choi, “Thermal Safety Analysis of a Dry Storage Cask for the 

Korean Standard Spent Fuel,” Paper No. ICEM2009-16159, Proceedings of the 12th

International Conference on Environmental Remediation and Radioactive Waste Management 

(ICEM 2009), Liverpool, United Kingdom, October 11–15, 2009, American Society of 

Mechanical Engineers, New York City, New York, 2009.
References

7-2

Coleman, H.W. and F. Stern, “Uncertainties in CFD Code Validation,” ASME Journal of Fluids 

Engineering, Vol. 119, pp. 795–803, 1997.

Das, K., D. Basu, J. Solis, and G. Zigh, “Computational Fluid Dynamics Modeling Approach to 

Evaluate VSC–17 Dry Storage Cask Thermal Designs,” Proceedings of the CFD for Nuclear 

Reactor Safety Applications (CFD4NRS-3) Workshop, Validation and Application of CFD and 

CMFD Codes to Nuclear Reactor Safety Issues, September 14–16, 2010, U.S. Nuclear 

Regulatory Commission, Washington, DC, 2010.

Elizade-Blancas, F., E. Karaismail, and I.B. Celik, “Application of Error Scaling Method in 

Conjunction with GCI,” Proceedings of the 3rd Workshop on CFD Uncertainty Analysis, 

October 23–24, 2008, Instituto Superior Técnico, Lisbon, Portugal, 2008.

Gomez, P., E. Araya, and M. Greiner, “Simulations of Heat Transfer Within the Fuel 

Assembly/Backfill Gas Region of Transport Packages,” PVP2005-71070, Proceedings of 2005 

ASME Pressure Vessels and Piping Division Conference, July 17–21, 2005, American Society 

of Mechanical Engineers, Denver, Colorado, 2005.

Harrington, P., “Design and Engineering Update,” Presentation to the U.S. Department of 

Energy, Nuclear Waste Technical Review Board, May 9, 2006, U.S. Department of Energy. 

McLean, Virginia, 2006.

Heng, X., G. Zuying, and Z. Zhiwei, “A Numerical Investigation of Natural Convection Heat 

Transfer in Horizontal Spent-Fuel Storage Cask,” Nuclear Engineering and Design, Vol. 213,

pp. 59–65, 2002.

Hoo, S.H., H.C. No, H.M. Kom, and E.H. Lee, “Full-Scope Simulation of a Dry Storage Cask 

Using Computational Fluid Dynamics,” Nuclear Engineering and Design, Vol. 240, 

pp. 4,111–4,122, 2010.

Kern, D.Q., Process Heat Transfer, McGraw Hill Publisher, New York City, New York, 1950.

McKinnon, M.A., R.E. Dodge, R.C. Schmitt, L.E. Eslinger, and G. Dineen, “Performance 

Testing and Analyses of the VSC–17 Ventilated Concrete Cask,” TR–100305, Electric Power 

Research Institute, Palo Alto, California, 1992.

Nishimura, M., H. Shibazaki, S. Fujii, and I. Maekawa, “Natural Convection Heat Transfer in the 

Horizontal Dry Storage System for the LWR Spent Fuel Assemblies,” Journal of Nuclear 

Science and Technology, Vol. 33, pp. 821–828, 1996.

NRC, “Identification and Prioritization of the Technical Information Needs Affecting Potential 

Regulation of Extended Storage and Transportation of Spent Nuclear Fuel,” ML120580143. 

U.S. Nuclear Regulatory Commission, Washington, DC, 2012a.

NRC, NUREG–2152, “Computational Fluid Dynamics Best Practice Guidelines for Dry Cask 

Applications,” U.S. Nuclear Regulatory Commission, Washington, DC, 2012b.

Oberkampf, W.L. and C.J. Roy, Verification and Validation in Scientific Computing. 

Cambridge University Press, Cambridge, Massachusetts, 2010.
References

7-3

Pelletier, D. and P.J. Roache, Chapter 13: Verification and Validation of Computational Heat 

Transfer, Handbook of Numerical Heat Transfer, 2nd Edition, W.J. Minkowycz, E.M. Sparrow, 

and Y. J. Murthy, editors, John Wiley & Sons, New York City, New York, 2006.

Perry, R.H. and D.W. Green, Perry’s Chemical Engineers’ Handbook, 8th Edition, McGraw-Hill 

Publisher, New York City, New York, 1997.

Phillips, T.S. and C. Roy, “Evaluation of Extrapolation-Based Discretization Error and 

Uncertainty Estimators,” AIAA-2011-215, Proceeding of the 49th AIAA Aerospace Sciences 

Meeting including the New Horizons Forum and Aerospace Exposition, American Institute of 

Aeronautics and Astronautics, Orlando, Florida, 4-7 January, 2011. 

Rezende, H.C., A.A. Santos, M.A. Navarro, and E. Jordao, “Verification and Validation of a 

Thermal Stratification Experiment CFD Simulation,” Nuclear Engineering and Design, Vol. 248,

pp. 72–81, 2012.

Richardson, L.F., “The Approximate Arithmatical Solution by Finite Differences of Physical 

Problems Involving Differential Equations, With an Application to the Stresses in a Masonry 

Dam,” Transactions of the Royal Society of London, Series A, Vol. 210, pp. 307–357, 1910.

Richardson, L.F. and J.A. Gaunt, “The Deferred Approach to the Limit,” Philosophical 

Transaction of the Royal Society of London, Series A, Vol. 226, pp. 299–361, 1927.

Roache, P.J., “Conservatism of the GCI in Finite Volume Computations on Steady State Fluid 

Flow and Heat Transfer,” ASME Journal of Fluids Engineering, Vol. 125, No. 4, pp. 731–732,

2003.

Roache, P.J., “Error Bars for CFD,” AIAA-2003-0408, AIAA 41st Aerospace Sciences Meeting 

and Exhibit, American Institute of Aeronautics and Astronautics, Reno, Nevada, 2003.

Roache, P.J. Verification and Validation in Computational Science and Engineering, Hermosa 

Publishers, Albuquerque, New Mexico, 1998.

Roache, P.J., “Quantification of Uncertainty in Computational Fluid Dynamics,” Annual Review 

of Fluid Mechanics, Vol. 29, pp. 123–160, 1997.

Roache, P.J., “Perspective: A Method for Uniform Reporting of Grid Refinement Studies,” 

ASME Journal of Fluids Engineering, Vol. 116, pp. 405–413, 1994.

Rohsenow, W.M., J.P. Hartnett, and Y.I. Cho, Handbook of Heat Transfer, 3rd Edition,

McGraw-Hill, New York, New York, 1998.

Shibazaki, H., M. Nishimura, N. Takahashi, S. Fujii, and I. Maekawa, “A Study of Heat Transfer 

Characteristics for a Horizontal Dry Storage System for LWR Spent Fuel Assemblies,” Heat 

Transfer-Japan Resources, Vol. 27, pp. 284–298, 1998.

Suffield, S.R., J.M Cuta, J.A. Fort, B.A. Collins, H.E. Adkins, and E.R. Siciliano, “Thermal 

Modeling of NUHOMS HSM-15 and HSM-1 Storage Modules at Calvert Cliffs Nuclear Power 

Station ISFSI,” PNNL–21788, Pacific Northwest National Laboratory, Richland, Washington,

October 2012.
References

7-4

TECPLOT, Inc, “TECPLOT 360 Version 2009,” TECPLOT, Inc., Bellevue, Washington, 

2009.Walavalkar, A.Y. and D.G. Schowalter, “3-D CFD Simulation of a Ventilated Concrete 

Cask Used for Spent Nuclear Fuel Storage,” Proceedings of the WM’04 Conference, February 

29 – March 4, 200, Tucson, AZ, Paper WM-4314.

Zigh, A. and J. Solis, “Computational Fluid Dynamics Best Practice Guidelines in the Analysis 

of Storage Dry Cask,” Proceedings of the Waste Management 2008 (WM2008) Conference, 

Phoenix, Arizona, February 24–28, 2008, WM Inc. Symposia, Tucson, Arizona, 2008.

Zoldners, N.G., “Thermal Properties of Concrete Under Sustained Elevated Temperatures,” 

ACI Publication, Paper SP25-1, American Concrete Institute, Detroit, Michigan, 1970.
APPENDIX A

NATURAL CONVECTION CORRELATION AND MODEL FOR VENT 

SCREENS AT INLET

A–1

APPENDIX A: NATURAL CONVECTION CORRELATION AND 

MODEL FOR VENT SCREENS AT INLET

A.1 Natural Convection Correlation

Convective heat transfer coefficients were specified at the exposed external surfaces of the 

horizontal storage modules HSM–1 and HSM–15. This appendix discusses how these 

coefficients were calculated using correlations for natural convection flows. All these 

correlations are expressed in terms of dimensionless parameters such as the Nusselt number 

(Nu), Rayleigh number (Ra), and Prandtl number (Pr).

The Nusselt number is the dimensionless ratio of the convection and conduction heat transfer 

coefficient and is defined as

Equation (A–1)

where

Nu = Nusselt number

h = Convection heat transfer coefficient

L = Length over which convection heat transfer takes place (characteristic length scale)

k = Thermal conductivity of the fluid.

The Prandtl number is a dimensionless number that is the ratio of momentum diffusivity 

(kinematic viscosity) to thermal diffusivity. It is expressed as

Equation (A–2)

where

Pr = Prandtl Number

ν = Kinematic viscosity of fluid

α = Thermal diffusivity of fluid.

The Rayleigh number is a dimensionless parameter that estimates the relative importance of 

convective heat transfer due to buoyancy-driven flow to heat conduction through the gas. It is 

expressed as 

Equation (A–3)

where

g = Gravitational acceleration

β = Thermal expansion coefficient

ΔT = Temperature difference between the heated surface and fluid

L = Characteristic Length.

The characteristic lengths used in the present study are the vertical or horizontal surface lengths

of the external surfaces of HSM–1 and HSM–15. This appendix provides the correlations to 

calculate the Nusselt number for both horizontal and vertical flat plates. These correlations are 

k

hL Nu =

α

ν Pr =

να

β∆ =

3 g TL Ra
Natural Convection Correlation and Model for Vent Screens at Inlet

A–2

obtained from Rohensow et al. [1998] and Bergman et al. [2011] and were derived for walls with 

uniform temperatures. Though the external wall surface temperatures of HSM-1 and HSM-15 

are not constant, the difference in the maximum and minimum temperature for a surface does

not exceed 5K [90 F]. Hence, application of these correlations provide a reasonable value of the 

heat transfer coefficient.

A.1.1 Correlations for Vertical Flat Plates

The Nusselt number for flow over a vertical flat plate with uniform wall temperature in laminar 

regimes is calculated as follows

ln[1 2.0 /(C Ra )]

2.0 Nu 1/ 4

l

l

+ = Equation (A–4)

where

Cl = 9 /16 4 / 9 [1 (0.492/ Pr) ]

0.671

+

Equation (A–5)

The Nusselt number in the turbulent region is given by

Equation (A–6)

where

Equation (A–7)

Once the laminar and turbulent Nusselt numbers are obtained, they are used to calculate the 

actual Nusselt number as 

[( ) ( ) ] m 1/ m

t

m Nu = Nul + Nu Equation (A–8)

where the constant m = 6 [Rohensow et al. 1998]. Equations (A–4) and (A–6) are valid for 10−1

< Ra < 1012.

A.1.2 Correlations for Horizontal Flat Plates

For horizontal plates with uniform wall temperature, the Nusselt number (Nu) for the laminar 

region is given by

ln[1 1.677 /(C Ra )]

1.4

ln[1 1.4/(0.835C Ra )]

1.4 Nu 1/ 4

l

1/ 4

l

l

+ = + = Equation (A–9)

is provided in Equation (A−5).

The Nusselt number for the turbulent regime is given by

Nu C (Ra) /(1 1.4 10 Pr/ Ra) 1/ 3 9

t = ν + ×

( )

0.42 0.81

0.22

1 0.61Pr

0.13Pr

+

Cν =

Cl
Natural Convection Correlation and Model for Vent Screens at Inlet

A–3

Equation (A–10)

Equations (A−9) and (A−10) are valid for 1 < Ra < 1010.

The actual Nusselt number is based on the laminar and turbulent Nusselt number and is 

given by

[( ) ( ) ] m 1/ m

t

m Nu = Nul + Nu Equation (A–11)

where the constant m = 10 [Rohensow et al. 1998].

A.2 Modeling of Vent Screens at Inlet

Both the HSM–1 and HSM–15 modules at the Calvert Cliffs Independent Spent Fuel Storage 

Installation (ISFSI) include a number of screens to protect the inlet and outlet vents. These 

screens will reduce the effective flow area and cause a pressure drop at the inlet and outlet. In 

the ANSYS-FLUENT [ANSYS 2012] models of HSM–1 and HSM–15, the inlets and outlets are 

modeled as vents, where a dimensionless pressure loss coefficient is specified to account for 

the resistance offered by the screens.

The pressure drop coefficient was calculated using the methodology Perry and Green [1997]

suggested for incompressible flows through screens. The following equation [Perry and Green 

1997] provides a functional relationship between pressure loss and momentum loss

Equation (A–12)

where

∆p = Pressure drop

ρ = Fluid density

V = Superficial velocity based upon the gross area of the screen

K = Pressure loss coefficient, given by

Equation (A–13)

where

C = Discharge coefficient

ε = Screen porosity.

The discharge coefficient is expressed as 

C = 0.1263 ln(Re) + 0.2043 (For Re<20) Equation (A–14)

C = 0.1 (Re)0.5 Equation (A–15)

Re is the screen Reynolds number 

( )

1/ 3

t Ra

1 0.01Pr

1 0.0107 Pr Nu 0.14 



 





+

+ =

2

V

p K

2 ρ ∆ =

















ε

− ε 



 



 = 2

2

2

1

C

1 K
Natural Convection Correlation and Model for Vent Screens at Inlet

A–4

Equation (A–16)

where 

Ds = aperture length 

μ = fluid viscosity

Equations (A–13) through (A–15) show that the pressure loss coefficient is a function of flow 

velocity. The pressure loss coefficient was calculated for a number of velocities, covering the 

range expected in the present study. In the calculations, it was assumed that the screens will 

reduce the flow area by 27 percent, and the aperture length of the screen is 0.0013716 m 

[0.054 in] [Suffield et al. 2012]. These pressure loss coefficients were specified as a piecewise 

linear function of normal velocities at the inlet and outlet vents of the ANSYS-FLUENT model.

εμ

ρ = D V Re s
APPENDIX B

MATERIAL PROPERTIES USED IN THE SIMULATIONS

B–1

APPENDIX B: MATERIAL PROPERTIES USED IN THE SIMULATIONS

Table B–1. Properties of air used in the simulations

[Bergman et al. 2011]

Temperature Specific heat Thermal conductivity Viscosity

J/kg-K BTU/lb-°F W/m-K Btu/[ft h °F] Pa-s lbf-s/ft2

250 K [−9.7 °F] 1,006 0.2402 0.0223 0.013 1.60 × 10−05 3.341 × 10−07

300 K [80.3 °F] 1,007 0.2405 0.0263 0.015 1.85 × 10−05 3.863 × 10−07

350 K [170.3 °F] 1,009 0.241 0.03 0.017 2.08 × 10−05 4.344 × 10−07

400 K [260.3 °F] 1,014 0.242 0.0338 0.019 2.30 × 10−05 4.803 × 10−07

450 K [350.3 °F] 1,021 0.244 0.0373 0.021 2.51 × 10−05 5.242 × 10−07

500 K [440.3 °F] 1,030 0.246 0.0407 0.023 2.70 × 10−05 5.639 × 10−07

600 K [620.3 °F] 1,051 0.251 0.0469 0.027 3.06 × 10−05 6.391 × 10−07

800 K [980.3 °F] 1,099 0.262 0.0573 0.033 3.70 × 10−05 7.727 × 10−07

Table B–2. Properties of helium used in the simulations

[Kern 1950]

Temperature Specific heat Thermal conductivity Viscosity

J/kg-K BTU/lb-°F W/m-K Btu/[ft h °F] Pa-s lbf-s/ft2

260 K [8.3 °F] 5,193 1.24 0.137 0.079 1.80 × 10−05 3.76 × 10−07

300 K [80.3 °F] 5,193 1.24 0.152 0.087 1.99 × 10−05 4.16 × 10−07

400 K [260.3 °F] 5,193 1.24 0.187 0.108 2.43 × 10−05 5.07 × 10−07

500 K [440.3 °F] 5,193 1.24 0.220 0.127 2.83 × 10−05 5.91 × 10−07

700 K [800.3 °F] 5,193 1.24 0.291 0.168 3.50 × 10−05 7.31 × 10−07

1,000 K [1,340.3 °F] 5,193 1.24 0.354 0.204 4.46 × 10−05 9.31 × 10−07

Table B–3. Properties of concrete used in the simulations

[Zoldners 1970]

Temperature Specific heat Thermal conductivity Density

J/kg-K BTU/lb-°F W/m-K Btu/[ft h °F] Kg/m3 lbm/ft3

311 K [100.13 °F] 1046.7 0.25 2.02 1.169 2,242 140

366.5 K [200.03 °F] 1046.7 0.25 1.97 1.139 2,242 140

533 K [499.73 °F] 1046.7 0.25 1.79 1.035 2,242 140

811 K [1,000.13 °F] 1046.7 0.25 1.38 0.798 2,242 140

Table B–4. Properties of carbon steel used in the simulations

[Kern 1950]

Temperature Specific heat Thermal conductivity Density

J/kg-K BTU/lb-°F W/m-K Btu/[ft h °F] Kg/m3 lbm/ft3

300 K [80.13 °F] 434 0.103 60.5 34.98 7,845 489.7

400 K [260.03 °F] 487 0.116 56.7 32.78 7,845 489.7

600 K [620.3 °F] 559 0.133 48 27.75 7,845 489.7

800 K [980.3 °F] 685 0.163 39.2 22.66 7,845 489.7

1,000 K [1,340.3 °F] 1,169 0.279 30 17.34 7,845 489.7
B–2

Table B–5. Properties of stainless steel used in the simulations

[Bergman et al. 2011]

Temperature Specific heat Thermal conductivity Density

J/kg-K BTU/lb-°F W/m-K Btu/[ft h °F] Kg/m3 lbm/ft3

300 K [80.13 °F] 480 0.114 15.1 8.73 8,055 502.8

400 K [260.03 °F] 512 0.122 17.3 10.00 8,055 502.8

600 K [620.3 °F] 559 0.133 20 11.56 8,055 502.8

800 K [980.3 °F] 585 0.139 22.8 13.18 8,055 502.8

1,000 K [1,340.3 °F] 606 0.144 25.4 14.68 8,055 502.8

Table B–6. Properties of lead used in the simulations

[Bergman et al. 2011]

Temperature Specific heat Thermal conductivity Density

J/kg-K BTU/lb-°F W/m-K Btu/[ft h °F] Kg/m3 lbm/ft3

300 K [80.13 °F] 129 0.0308 35.3 20.41 11,340 707.9

400 K [260.03 °F] 132 0.0315 34 19.65 11,340 707.9

600 K [620.3 °F] 142 0.0339 31.4 18.15 11,340 707.9

Table B–7. Emissivity of materials used in the simulations

[Bergman et al. 2011]

Material Emissivity

Stainless Steel 0.46

Concrete 0.9

Carbon Steel 0.65
APPENDIX C

POROUS MEDIA MODEL

C–1

APPENDIX C: POROUS MEDIA MODEL

Dry storage canister (DSC) shells hold 24 pressurized water reactor (PWR) fuel assemblies in 

both the HSM–1 and HSM–15 horizontal storage modules. Each of these fuel assemblies 

houses an arrangement of Combustion Engineering (CE) 14 × 14 fuel rods. The most accurate 

method for simulating the fuel assembly would involve explicitly modeling each fuel rod in the 

assembly, along with all the accessories such as spacer grids and guide thimbles. However, 

such explicit modeling would be computationally intensive. A widely employed and realistic 

alternative is to use a continuous porous medium to represent the volume within a fuel 

assembly instead of modeling all the internal components. In general, the helium flow velocity is 

expected to be low within the fuel assemblies and characterized by a large pressure drop

across the length of the assembly due to the presence of a large number of solid obstacles. In 

the porous media approximation, predetermined flow resistances are specified in each flow 

direction to simulate the effect of pressure loss due to the presence of the fuel rods and other 

obstacles. 

In ANSYS-FLUENT, the porous media model uses a volume-based source term in the 

momentum equation to account for the pressure loss throughout the porous zone [ANSYS Inc. 

2012]. The pressure drop due to flow resistance in the porous region is expressed in the 

following equation

Equation (C–1)

where

μ = Viscosity of helium

ρ = Density of helium

V = Velocity magnitude in axial direction

p = Pressure

The first term in Equation (C–1) represents viscosity-dominated (frictional) losses and the 

second term provides the inertia-dominated loss. The constants C and D are known as inertial 

and frictional flow resistance factors, respectively. For the present study, these flow resistance 

factors need to be specified for the actual configuration of the fuel assembly. 

To obtain the inertial and frictional flow resistance parameters C and D, a separate analysis was 

performed using a three-dimensional (3-D) computational fluid dynamics (CFD) simulation of an 

isolated fuel assembly. The computational domain consisted of a quadrant of the assembly to 

take advantage of the symmetry in the fuel arrangement. In the axial direction, the 

computational domain spanned the entire assembly length. Figure C–1 shows a threedimensional (3-D) view of the entire computational domain. Fuel and control rods were explicitly 

modeled as solid cylindrical objects. Figure C–2 provides a closer look into the domain and 

shows the arrangement of fuel and control rods in the quadrant. Table C–1 lists the pertinent 

dimensions used in this analysis. The computational domain was meshed using a hexahedral 

mesh, as shown in Figure C–3. At the inlet, a velocity boundary condition was used. At the 

outlet, a pressure of one atmosphere—that is, a zero gauge pressure—was specified.

2 V

2

1 D V C

dx

dp = • μ + • ρ  
Porous Media Model

C–2

Figure C–1. Longitudinal full view of the computational domain for the CE 

14 × 14 fuel assembly for porous media parameter calculation showing the 

fuel rod

Figure C–2. 1/4th cross-sectional domain for the CE 14 × 14 fuel assembly for porous 

media parameter calculation
Porous Media Model

C–3

Table C–1. Dimensions used in the simulations

Length of the tube 4.1275 m [13.54 ft]

Fuel rod diameter 0.011176 m [0.037 ft]

Guide tube diameter 0.028321 m [0.093 ft]

Cross sectional area at inlet/outlet 0.007695 m2 [0.083 ft2]

Figure C–3. Two-dimensional view of the computational domain for the CE 

14 × 14 fuel assembly for porous media parameter calculation

From Equation (C–1), it can be seen that the pressure drop is a function of velocity through the 

porous medium. Hence, to determine the frictional and inertial loss factors, the simulations 

need to be performed for at least two inlet velocity sets. Based on the computed pressure drops 

obtained from these simulations, the loss factors can be determined. For the present study, a 

series of inlet velocities was used to cover the range expected under normal operating 

conditions for the storage system. Using several pairs of the inlet velocity–pressure drop sets

listed in Table C–2 , a set of values for C and D in equation C–2 was determined. The average 

values of the calculated C and D from this set were used in the simulations. The calculated 

average inertial (C) and viscous pressure loss factors (D) are 0.78 1/m (0.2377 1/ft) and 

1.2848 × 105 1/m2 (3730.2883 1/ft2

), respectively. Note that these factors are only for the axial 

flow direction. Due to the presence of the wall and a large number of obstructions, it was 

assumed that the flow resistances in the transverse directions are much higher. Consequently, 

arbitrarily high inertial and viscous pressure loss factors (two orders of magnitude higher than

the axial direction) was used for the transverse directions. A sensitivity study, reported in 

Chapter 3, was conducted to assess the effect of the resistance factors.
Porous Media Model

C–4

Table C–2. Inlet condition and parameters for calculation of loss factors

Inlet velocity Pressure drop (∆p)

0.01 m/s [0.032 ft/s] 0.1056 Pa [2.206 × 10−3 lbf/ft2]

0.02 m/s [0.064 ft/s] 0.2112 Pa [4.411 × 10−3 lbf/ft2]

0.04 m/s [0.128 ft/s] 0.4226 Pa [8.826 × 10−3 lbf/ft2]

0.06 m/s [0.192 ft/s] 0.6342 Pa [1.325 × 10−2 lbf/ft2]

0.08 m/s [0.260 ft/s] 0.8461 Pa [1.767 × 10−2 lbf/ft2]

0.10 m/s [0.32 ft/s] 1.0582 Pa [2.21 × 10−2 lbf/ft2]
APPENDIX D

CALCULATION OF FUEL EFFECTIVE THERMAL CONDUCTIVITY

D–1

APPENDIX D: CALCULATION OF FUEL EFFECTIVE THERMAL 

CONDUCTIVITY

The fuel assemblies inside the dry storage canister (DSC) shell contain a complex arrangement 

of fuel rods and associated accessories. It would be computationally intensive to explicitly 

model every fuel rod; therefore, the fuel assemblies are modeled as a continuous porous zone 

having an anisotropic thermal conductivity. The effective thermal conductivity will represent the 

effective conduction and radiation heat transfer in the fuel assemblies. The anisotropy in 

thermal properties is caused by significant differences in thermal conductivities in radial and 

axial directions. The thermal conductivity in the radial direction is interrupted by gaps between 

rods that are occupied by helium, whereas in the axial direction, heat flow primarily occurs 

through the fuel cladding walls. Heat transfer between the fuel rods and assembly walls is 

dominated by radiation heat transfer, which is dependent on the temperature difference 

between participating walls, so the effective thermal conductivity is a function of temperature. 

This appendix describes the analysis that was done to calculate the axial and effective thermal 

conductivities for the Combustion Engineering (CE) 14 × 14 fuel rods under the storage 

conditions expected in HSM-1 and HSM-15. 

D.1 Radial Thermal Conductivity 

To calculate the radial effective thermal conductivity, a two-dimensional (2-D) ANSYS-FLUENT 

[ANSYS 2012] model was developed to represent an axial cross section of the fuel assembly 

with fuel rods. The 2-D model is based on a quarter cross section of the CE 14 × 14 assembly 

and includes the full geometrical details of the fuel rods and Zircaloy cladding. The fuel rods 

were modeled as a homogeneous solid material region with a specified, uniform heat generation 

rate. Thermal conductivities of the UO2 fuel and Zircaloy were used for the fuel and cladding 

region, respectively. The interstitial space between the fuel rods and the assembly walls was

assumed to be occupied by the helium backfill gas and, consequently, the thermal conductivity 

of helium is used in the interstitial space. The computational domain and the grids used in the 

2-D simulations are shown in Figure D–1. A uniform, volumetric heat generation rate {8,522 

W/m3 [824.03 BTU/hr ft3

]} was applied to the fuel region, and uniform temperature was applied 

to the edges of the 2-D computational domain. The volumetric heat generation rate used in the 

current analysis was calculated by dividing the total decay heat load of the fuel region by the 

equivalent volume of the fuel region. Radiation was modeled using the Discrete Ordinate (DO) 

model. Symmetry conditions were used for the left and lower sides of the domain. 

Temperature values were specified for the top and right sides of the domain, which represent 

the outer bound of the fuel assembly. 

A number of steady-state simulations were performed for a range of boundary temperatures: 

300, 350, 400, and 900 K [80.3, 170.3, 260.3, and 1,160.3 °F]. Figure D–2 shows the 

temperature distribution for a simulation with an edge temperature of 300 K [80.3 °F]. Based on 

the temperature distribution and temperature gradient obtained from the series of steady-state 

simulations, the fuel effective thermal conductivity in the radial direction was determined using 

the keff methodology described by Bahney and Lotz [1996]. The values of effective radial 

thermal conductivities for the fuel region are listed in Table D–1. These values were used in the 

ANSYS-FLUENT model as an orthotropic temperature-dependent radial conductivity within the 

fuel assembly, approximated as a porous medium.
Calculation of Fuel Effective Thermal Conductivity

D–2

Figure D-1. Computational domain and grid for fuel, clad, and helium gas grids in the 

effective thermal conductivity calculation

Figure D–2. Predicted temperature contours for an edge temperature of 300 K [80.3 °F]
Calculation of Fuel Effective Thermal Conductivity

D–3

Table D–1. Keff of fuel region in the radial direction with helium backfill

Temperature (K) [°F] Keff (radial) (W/m-k) [Btu/(ft h °F)]

307.09 K [93.1 °F] 0.3728 [0.22]

405.08 K [269.4 °F] 0.5203 [0.30]

503.70 K [447 °F] 0.7144 [0.41]

602.74 K [625.26 °F] 0.9647 [0.56]

D.2 Axial Thermal Conductivity 

The effective thermal conductivity in the axial direction was determined using an area-weighted 

average of the thermal conductivities of the materials within the fuel area (i.e., KUO2, Khelium, and 

KZircaloy). The relationship used to calculate the effective thermal conductivity in the axial 

direction is shown in Equation (D–1).

Keffective = KUO2 × Aratio−UO2 + Khelium × Aratio−helium + KZircaloy × Aratio−Zircaloy Equation (D–1)

Equation (D–1) uses the area ratio of each material (i.e., Aratio-UO2, Aratio-helium, and Aratio-Zircaloy) with 

respect to the total area of the cross section of the fuel region and calculates the effective 

thermal conductivity in the axial direction. The calculated values of temperature dependent 

effective axial thermal conductivity are presented in Table D–2.

These listed values were specified as temperature-dependent orthotropic effective 

conductivity in the axial direction within the fuel assembly volumes, approximated as 

continuous porous media.

Table D–2. Keff of fuel region in the axial direction with helium backfill

Temperature (K) [°F] Keff (axial) (W/m-k) {Btu/[ft h °F]}

300 K [80.33 °F] 2.5689 [1.49]

400 K [260.33 °F] 2.6250 [1.52]

500 K [440.33 °F] 2.7141 [1.57]

600 K [620.33 °F] 2.7656 [1.6]]

APPENDIX E

FUEL ASSEMBLY LOADING

E–1

APPENDIX E: FUEL ASSEMBLY LOADING

The decay heat load was defined for each of the 24 active fuel assemblies within the dry 

storage canister (DSC) in horizontal storage modules HSM–1 and HSM–15 based on decay 

heat load data from the Calvert Cliffs Independent Spent Fuel Storage Installation (ISFSI). The 

fuel structure and fuel identification location for HSM–1 and HSM–15 are shown in Figures E–1 

and E–2, respectively. The total decay heat load for the DSC in the HSM–1 module is 

calculated at approximately 4.1 kW [13,990 BTU/hr] as of June 2012. The total decay heat load 

for HSM–15 as of June 2012 is 7.6 kW [25,933 BTU/hr]. Tables E–1 and E–2 show the decay 

heat load values for each fuel region in the HSM–1 and HSM–15 modules as of June 2012. 

This value of the heat load is divided by the volume of the active fuel region for each assembly 

and used in the simulations. 

The decay heat load was applied as a uniform volumetric heat generation rate over the active 

fuel length within each fuel region. The volumetric heat generation rate was modified to include 

the axial variation. This was done by multiplying the volumetric heat generation rate by an axial 

peaking factor. The peaking factor was based on Combustion Engineering (CE) 14 × 14 

representative burn up data from the Calvert Cliffs ISFSI [Suffield et al. 2012]. The data were

obtained from a 2006 study at Calvert Cliffs to determine the thermal and radiological source 

terms for the standard CE 14 × 14 fuel assemblies in the spent fuel pool. The axial variation of 

peaking factor was implemented as a combination of linear segments in the energy source term 

specified in the ANSYS-FLUENT model [ANSYS Inc. 2012]. Figure E–3 shows the axial 

distribution of the peaking factor against nondimensional axial length (length of fuel region). 

This distribution is used in all the simulations for HSM–1 and HSM–15. 
Fuel Assembly Loading

E–2

Figure E–1. Fuel assembly loading identification structure for HSM–1 module

Figure E–2. Fuel assembly loading identification structure for HSM–15 module
Fuel Assembly Loading

E–3

Table E–1. Calculated decay heat load values at year 2012 for HSM–1 module

Assembly ID (fuel region)

for HSM-1 module

2012 heat load

(Current year of measurement)

(Watt)

1 Watt = [3.413 BTU/hr]

1A007 130.5

1A034 130.7

1A035 133.0

1A023 137.5

1A041 137.4

1A059 137.4

1A029 138.2

1A054 139.3

1B032 139.9

1B013 183.9

1B053 184.7

1B079 187.4

1B046 192.2

1B044 195.5

1B004 197.8

1B018 197.7

1B061 198

1B003 198.9

1C005 175.2

1C103 189.9

1C107 189.4

1C106 198.3

1C206 202.1

1C214 201.7

Table E–2. Calculated decay heat load values at year 2012 for HSM–15 module

Assembly ID

for HSM–15 module

2012 heat load

(current year of measurement)

(Watt)

1 Watt = 3.413 BTU/hr

2F177 272

2F136 291

1E 101 301

2F021 385

2F024 385

2F118 286

2F179 290

1H115 406

2G123 427

1H109 407

1G018 327

2F117 285

2F123 285

1G027 340
Fuel Assembly Loading

E–4

Table E–2. Calculated decay heat load values at year 2012 for HSM–15 module (continued)

Assembly ID

for HSM–15 module

2012 heat load

(current year of measurement)

(Watt)

1 Watt = 3.413 BTU/hr

1H119 408

1H124 407

1E 102 301

1F119 232

2F166 272

1E 121 302

1E 120 302

1F110 225

2F125 274

1F113 225

Figure E–3. Axial distribution of the peaking factor used in current simulations





UNITED STATES

NUCLEAR REGULATORY COMMISSION

WASHINGTON, DC 20555-0001

------------------

OFFICIAL BUSINESS
NUREG/CR

-7191 Thermal Analysis of Horizontal Storage Casks for 

Extended Storage Applications