NUREG-1726

Predictions of Spent Fuel 

Heatup After a Complete 

Loss of Spent Fuel Pool 

Coolant 

U.S. Nuclear Regulatory Commission O 

Office Nuclear Regulatory Research 

Washington, DC 20555-0001
AVAILABILITY OF REFERENCE MATERIALS 

IN NRC PUBLICATIONS

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 www.nrc.gov/NRC/ADAMS/index.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 

P. O. Box 37082 

Washington, DC 20402-9328 

www.access.gpo.gov/sudocs 

202-512-1800 

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: Office of the Chief Information Officer, 

Reproduction and Distribution 

Services Section 

U.S. Nuclear Regulatory Commission 

Washington, DC 20555-0001 

E-mail: DISTRIBUTION@nrc.gov 

Facsimile: 301-415-2289 

Some publications in the NUREG series that are 

posted at NRC's Web site address 

www.nrc.gov/NRC/NUREGS/indexnum.html 

are updated regularly and may differ from the last 

printed version.

Non-NRC Reference Material 

Documents available from public and special 

technical libraries include all open literature 

items, such as books, journal articles, and 

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 4 2 nd Street 

New York, NY 10036-8002 

www.ansi.org 

212-642-4900

The NUREG series comprises (1) technical 

and administrative reports and books prepared 

by the staff (NUREG-)QQQ() or agency 

.contractors (NUREG/CR-XXXX), (2) ...proceedings of conferences 

(NREG/CP-XXXX), (3) reports resulting from 

:international agreements (NUREG/IA-XXXX), 

(4) brochures (NUREGIBR-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);
NUREG-1726

Predictions of Spent Fuel 

Heatup After a Complete 

Loss of Spent Fuel Pool 

Coolant

Manuscript Completed: June 2000 

Date Published: June 2000 

Prepared by 

C.F. Boyd 

Division of Systems Analysis and Regulatory Effectiveness 

Office of Nuclear Regulatory Research 

U.S. Nuclear Regulatory Commission 

Washington, DC 20555-0001
ABSTRACT

This generic study was done to predict the peak fuel temperatures and flow patterns observed 

after a complete loss of spent fuel pool coolant. The pool is assumed to be housed in a large 

containment building with an operable ventilation path. Fuel temperatures and natural 

circulation phenomena are predicted using computational fluid dynamics (CFD). These 

predictions identify the peak temperatures and the principal flow paths and features of the 

natural convective cooling phenomena for the coupled system of the fuel pool and containment 

building. The model is sized to represent a typical BWR pool and containment building. The 

pool is filled to its capacity with 4200 fuel bundles of various ages in high-density racking. 

Predictions are made with fuel loads representative of a fuel pool 2, 3, 4, and 6 years after the 

reactor is shut down. Sensitivity studies are done on the ventilation rate, the outer wall heat 

transfer coefficient, the location of the hottest fuel, fuel burnup, the flow resistance within the 

racks, and heat conduction within the racks. Best estimate predictions of critical decay time 

show that fuel temperatures remain below the temperature limits of 800 °C and 600 °C after 26 

and 35 months of decay time, respectively. Predictions are made using the FLUENT CFD 

code. The finite volume mesh used for this model provides the resolution needed to resolve the 

important phenomena while remaining small enough to provide solutions with the available 

computer equipment. The overall size and computational complexity of the completed model 

preclude a rigid grid independence study. In considering quantitative results, the limitations and 

modeling assumptions must be kept in mind.

i°i
CONTENTS

Page 

ABSTRACT iii 

FIGURES vi 

TABLES vii 

1 INTRODUCTION 1 

2 COMPUTATIONAL FLUID DYNAMICS (CFD) 2 

3 PHYSICAL MODEL DESCRIPTION 2 

3.1 Pool and Racks 2 

3.2 Fuel 3 

3.3 Containment Building 4 

3.4 Ventilation 5 

3.5 Wall Heat Transfer 5 

4 CFD MODEL DEVELOPMENT 5 

4.1 Finite Volume Mesh 6 

4.2 Boundary Conditions 6 

4.3 Material Properties 7 

4.4 Fuel and Racks 8 

4.5 Turbulence Modeling 8 

4.6 Solution Convergence 9 

4.7 Grid Independence 9 

5 COMPLETED PREDICTIONS 9 

6 SUMMARY OF ASSUMPTIONS AND LIMITATIONS 11 

7 RESULTS 12 

7.1 Base Case Results (f4V2hO) 12 

7.2 Tabulated Results for Each Case 21 

7.3 Decay Time Sensitivity, Adiabatic Walls 23 

7.4 Ventilation Rate Sensitivity 25 

7.5 Wall Heat Loss Sensitivity 26 

7.6 Rack and Fuel Flow Resistance Sensitivity 27 

7.7 Hottest Fuel Location Sensitivity 28 

7.8 Rack and Fuel Thermal Conductivity Sensitivity 29 

7.9 Fuel Burnup Sensitivity 29 

7.10 Best-Estimate Critical Decay Time 30 

8 SUMMARY 30 

9 REFERENCES 32 

A APPENDIX A LOSS COEFFICIENTS FOR FLOW THROUGH A-1 

BWR FUEL STORAGE RACKS 

B APPENDIX B THERMAL PROPERTY DETERMINATION FOR POROUS B-1 

REGION REPRESENTING BWR FUEL AND STORAGE RACKS

V
FIGURES 

1. Pool Region and Racks 3 

2. High-Density Fuel Rack 3 

3. Containment Building With Fuel Pool 4 

4. Finite Volume Mesh for CFD Predictions 6 

5. Temperature Contours on Active Fuel Region 12 

6. Temperature Contours Showing Rising Hot Plume 13 

7. Temperature Contours Showing Cooling Air Flow 14 

8. Velocity Vectors in a Plane 0.5 m Above the Containment Floor 14 

9. Vertical Velocity (mis), Horizontal Plane, Top of Active Fuel 15 

10. Velocity Vectors Shaded With Temperature Below the Fuel Racks 16 

11. Pressure Contours Above the Fuel Racks, Pa (z = 4.5 m) 18 

12. Pressure Contours Below the Fuel Racks, Pa (z = 0.1 m) 19 

13. Temperature Contours Above the Fuel Racks (z = 4.5 m) 19 

14. Temperature Contours Below the Fuel Racks (z = 0.1 m) 20 

15. Predicted Temperatures vs. Decay Time, Adiabatic Walls 24 

16. Sensitivity of Temperatures to the Ventilation Rate 25 

17. Sensitivity of Temperatures to Overall Heat Transfer Coefficient 26 

18. Sensitivity of Temperatures to Flow Resistance 28 

19. Best-Estimate Results for Maximum Temperature vs. Decay Time 30

vi
TABLES 

Page 

1. Fuel Loads for Whole Pool 4 

2. Material Properties, Air (White, 1988) 7 

3. Completed Predictions 10 

4. Assumptions and Limitations 11 

5. Mass Flow Rates at Top of Active Fuel, Base Case 17 

6. Predicted Temperatures and Pressures 21 

7. Predicted Mass Flow Rates 22 

8. Predicted Global Energy Balance Terms 23

vii
1 INTRODUCTION

In support of the Nuclear Regulatory Commission (NRC) rulemaking activity related to 

decommissioning, the office of Nuclear Reactor Regulation (NRR) is completing a study on 

spent fuel pool accident risks. In support, the Office of Nuclear Regulatory Research (RES) is 

providing technical assistance in several areas. This report documents a thermal-hydraulic 

evaluation of spent fuel pool heatup after a low-probability complete loss of spent fuel coolant. 

Computational fluid dynamics (CFD) is used to predict fuel heatup and natural circulation flow 

paths throughout the spent fuel pool and the upper containment building. The predictions give 

insights into the phenomena of air cooling which provides most of the heat removal capacity 

during long-term cooling scenarios after a complete loss of coolant. 

Spent fuel pool heatup predictions are typically made using codes tailored to the geometry and 

physics of spent fuel stored in a rack. Codes of this type include SHARP 1, SFUEL, and 

COBRA-SFS. The SHARP manual (Ref. 1) provides general background information on some 

of the important phenomena considered in this type of model. A typical modeling approach is to 

determine one-dimensional buoyancy-driven flows between idealized upper and lower control 

volumes tying the fuel bundles together. Codes such as COBRA-SFS and SFUEL incorporate 

physical models for conduction, radiation, flow losses, clad oxidation chemistry, and other 

things. The flowfield assumptions, however, simplify the natural circulation flow paths in and 

around the fuel racks. The assumption of a single well-mixed volume joining each of the 

bundles at the top and bottom of the racks does not account for pressure or temperature 

variations in this region. The idealized upper control volume provides the ultimate heat sink for 

these models, simulating the entire building and the ventilation system. 

Previous studies indicate that in steady-state conditions, the heat produced by the fuel is 

removed primarily through natural convection flows 2 (Ref 2, p. 47). In addition, the largest 

source of uncertainty in these fuel heatup predictions is the natural circulation flow rate 3 (Ref 3, 

p. 57). The RES predictions highlight the natural circulation flows that are simplified in typical 

spent fuel pool models. The focus of the predictions is on the three-dimensional natural 

circulation flowfield in and around the fuel pool, racks, and containment building. Physical 

models for radiation and clad chemistry are not incorporated. The predictions can be used to 

assess the flowfield assumptions used in other codes. The three-dimensional CFD predictions 

made by RES give valuable insights into the natural circulation air flow which is crucial to spent 

fuel pool cooling after a complete loss of spent fuel pool coolant. 

All code results must be viewed in the light of the assumptions used for the analysis. The CFD 

model makes simplifying assumptions to represent the complex geometry of the racks and fuel. 

A porous region provides an equivalent flow resistance for the rack and bundles and aligns the 

flow in the vertical direction. A volumetric heat source adds the appropriate energy to the fluid 

in the active fuel region. Predictions are obtained for steady-state conditions to determine the 

maximum fuel surface temperature for a given pool age and configuration. The CFD 

predictions do not include models for radiation and clad chemistry. These models are 

considered important at elevated temperatures (T > 600 'C, Ref 2, Figs. 12, 14 and 26). 

Therefore, the current CFD predictions are more applicable at low temperatures, where these 

effects are minimal. At elevated temperatures, potential limitations should be kept in mind.

I
2 COMPUTATIONAL FLUID DYNAMICS (CFD)

The FLUENT4 (version 5) CFD code is used to predict the temperature rise in the fuel and the 

natural circulation flow rates. FLUENT is a commercially available, general-purpose CFD code 

capable of solving a wide variety of fluid flow and heat transfer problems. The code solves the 

Reynolds-averaged Navier-Stokes equations a finite volume mesh. The Navier-Stokes 

equations are the mass, momentum, and energy equations for a continuous fluid. Reynolds

averaging creates the need for turbulence modeling to account for the turbulent diffusion of 

momentum and energy. The FLUENT code provides several turbulence modeling options and 

can be applied to almost any geometry. 

The predicted results from a CFD analysis are influenced by several steps in the process. In a 

blind study, where no experimental data are used for bench marking the code, the analyst must 

rely on accepted CFD practices and engineering judgment to create an acceptable CFD 

simulation. The basic steps of a CFD analysis are to describe the physical model, develop the 

CFD model, do sensitivity studies, and validate the solution. Commercial CFD codes are widely 

used in many industries today. The fundamentals of CFD are to be found in textbooks such as 

the introductory text by Anderson5 . 

3 PHYSICAL MODEL DESCRIPTION 

A simplified model of the containment building, fuel pool, racks, and ventilation system is 

developed for this investigation. The major parameters used and the physical features of the 

pool are determined with input from the Office of Nuclear Reactor Regulation (NRR). The 

model is based upon a boiling water reactor (BWR) containment and fuel pool. The major 

components are described below. 

3.1 Pool and Racks 

The pool and fuel rack dimensions are typical of a large BWR. The pool is assumed to be 

completely filled with high-density racking. The racks extend to within 0.15 m (6 in.) of the fuel 

pool Walls. To establish symmetry, the pool is centered on the long axis of the containment 

building. Figure 1 shows the features of the fuel pool cut along the line of symmetry. The finite 

volume cells used to represent the rack region in the CFD model are included to highlight the 

location of the racks. Overall pool dimensions are 9.54 m x 11.08 m x 12 m (I, w, d). The fuel 

racks sit 0.15 m (5.9 in.) from the pool floor and are 4.3 m (14.1 ft) in height. The high-density 

racking occupies the entire lower pool (except the 0.15 m gap around the exterior). Figure 2 

shows the structure of a typical high-density rack. The square cells which make up the rack 

share a wall between neighboring cells and there are no open cells. The walls are solid so 

there is no communication of fluid between cells. Cell pitch is 0.154 m (6 in.). Rack-to-rack 

spacing is ignored in this model, which is equivalent to assuming a single rack structure with 

4200 cells (60 x 70; 9.24 m x 10.78 m). Assuming a full core contains 800 bundles, this pool is 

racked to hold 5 1/4 cores.

2
'symmetry plane 

finite volumes from 

CFD model shown 

for racks

Figure 1. Pool Region and Racks Figure 2. High-Density Fuel Rack

A fuel bundle is assumed to occupy each cell of the high-density racking (4200 assemblies). A 

uniform active (heat generating) length of 3.8 m (12.46 ft) is centered from top to bottom in the 

racks. Four pool regions are defined with unique characteristics. A hot region with 800 bundles 

(fl = 1 core), a second region with 267 bundles (f2 = 1/3 core), a third region with 267 bundles 

(f3), and a final region with 2866 bundles (f4). Fuel f2 is assumed to be 18 months older than 

fuel fl and fuel f3 is 18 months older than fuel f2. Region f4 represents all of the remaining 

(older) fuel. The fuel layout is illustrated in Figure 1. The hot fuel (fl) is assumed to remain in 

a contiguous region on the left side of the pool. Fuel f2 is in another contiguous region next to 

the hottest fuel. Next comes fuel f3. The remaining 2866 bundles, f4, fill out the rack. A 

sensitivity study is done to look at the effect of moving the hottest fuel, fl, to the center of the 

pool. More details on the geometry of the 9 x 9 bundles are given in Appendix A, "Flow Loss 

Coefficients through BWR Fuel Storage Racks." Table 1 lists the fuel loads applied in this 

analysis. The fuel loads are given in watts per bundle. These loads are based upon an 

average burnup of 40 or 50 GWd/MTU, as noted. Most of the predictions assume an average 

burnup of 40 GWd/MTU for each fuel region. A sensitivity study is done assuming a burnup of 

50 GWd/MTU.

3

3.2 Fuel
Table 1. Fuel Loads for Whole Pool

No. of No. of Average Watts/bundle 

cores cells after number of years 

2 yr 3 yr 4 yr 4yr 6 yr (40)* (40)* (40)* (50)* (40)* 

fl, Hot 1 800 920 620 462 576 341 

f2, Med Hot .33 267 536 419 360 457 296 

f3, Med Cold .33 267 381 341 304 390 272 

f4, Cold 3.6 2866 201 198 194 246 190 

Total heat load (megawatts) >>> 1.557 1.266 1.103 1.392 0.969 

* Average Burnup (GWd/MTU) 

3.3 Containment Building 

The main building is sized to represent a BWR containment. No internal structures are 

modeled. The size and location of the containment relative to the fuel pool is shown in Figure 

3. This figure represents the CFD model domain, which is a symmetrical half of the physical 

structure. The plane of symmetry faces out toward the reader. The containment building is 

modeled only from the pool surface up to the ceiling. Lower levels are ignored.

•symmetry plane

Main Building 

(1/2 shown)

2i-nMto pool edge ~--

2m x 2m output 

vent located on end wall

pool region

rack region 

Figure3.Containment Wih...... Building. ..Fl P 

Figure 3. Containment Building With Fuel Pool

4

ventilation 

input 2m x 2m
3.4 Ventilation

A single ventilation input and output are used in the CFD model, which represents half the 

building and fuel pool. The ventilation opening is 4 m2 and is shown in Figure 3. The standard 

ventilation rate is based upon 2 building volumes of air exchanged per hour (Ref. 2, p. 41). The 

inlet flow provides 2 building volumes of air to the building every hour. Sensitivity studies are 

performed on this value. Air enters the containment through the inlet vent and is heated by the 

fuel. Inlet air is set at 27 0C. Air is exhausted on the opposite end of the building through a 

similar vent opened to the atmosphere. The inlet and outlet vents are positioned 2.68 m from 

the back wall and ceiling of the containment building. 

3.5 Wall Heat Transfer 

Most of the predictions assume the external building walls are adiabatic. Sensitivity studies are 

done on the overall heat transfer coefficient at the side walls and ceiling of the containment 

building. Values of 0, 2, and 4 W/m 2-K are used. Best estimate values range6 from 1 to 3 

W/m2-K. This range covers various wall types subjected to a 15 mph outdoor wind speed. 

Outdoor air is assumed to be 27 0C. A value of 2 W/m 2-K is used for best estimate predictions 

in this report. The pool walls and containment floor are always modeled as adiabatic. 

4 CFD MODEL DEVELOPMENT 

The physical model described above is represented on a set of finite volumes on which the 

governing equations are discretized and solved. The Reynolds-averaged Navier-Stokes 

equations are used. Boundary conditions are applied on the model to facilitate a solution. 

Developing a CFD model requires a balance between the size and complexity of the model and 

the need to obtain solutions within available time and computational resources. Although the 

computational resources available to RES are significant, the scale of this problem dictates that 

simplifications must be made. The assumptions used to develop the CFD model are discussed 

below. The model is for use in the FLUENT version 5 unstructured CFD code. 

Each prediction is run to find a steady-state solution, using typical CFD modeling options. The 

solutions give the maximum expected temperature but not the time it takes to reach it. A 

transient analysis would require significantly more computer time. The major features of the 

CFD model are given below. 

steady state Reynolds averaged Navier-Stokes solution 

k-epsilon turbulence model with standard wall functions 

ideal gas, variable thermal properties (air as working fluid) 

gravity 

675,000+ finite volume cells (tetrahedron and hexagonal mesh) 

segregated solver with 2nd order differencing on momentum and energy 

porous media model of fuel and racks (vertical laminar flow through region) 

The clad oxidation reaction and radiative heat transfer are not included in this model. The 

effect of these omissions is considered small at temperatures below 600 0C. Higher 

temperature results should be considered with this limitation in mind. These effects may be 

studied in the future using the FLUENT code.

5
4.1 Finite Volume Mesh 

The mesh consists of over 675,000 finite volumes (Figure 4). The mesh design minimizes the 

number of cells used for the model while providing an adequate number of cells to resolve 

important physical gradients. Cells are constructed to have minimal skew and sides of equal 

length (aspect ratio = 1) wherever possible. Cell to cell growth is limited to a range between 0 

% and 20 %. 

Cells in the upper building are large, 1 M3 , to save resources, so excess diffusion is expected 

in this region. This concession is not considered serious. Flows in the upper building are 

expected to be more complex (more diffusive) in the typical containment due to internal 

structures omitted in the CFD model. Any numerical diffusion resulting from the large cells in 

the upper building will enhance mixing, which is expected to be underpredicted by this model.

'fl *22

Figure 4. Finite Volume Mesh for CFD Predictions 

Most of the 675,000 cells are located in and around the fuel racks where the highest gradients 

and most significant phenomena occur. Cell height at the top of the racks is reduced to 1 cm 

through the use of grid adaption. This fine resolution is needed to resolve the large thermal and 

velocity gradients predicted in this region. Transitions from the very fine mesh to regions with 

larger mesh are constructed of tetrahedral elements with a controlled growth rate. This 

technique eliminates the skew in the hexagonal element regions and provides a smooth 

transition between cells of different sizes. 

4.2 Boundary Conditions 

Boundary conditions are applied at the walls (both internal and external) and at the inlets and 

outlets. All external walls use the no-slip condition for the momentum equations. No surface 

roughness is applied. All external walls are adiabatic except as noted in the sensitivity studies.

6

•!•.
The surface facing the reader in Figure 4 is a symmetry plane. Internal vertical walls are 

positioned within the fuel and rack region to separate zones of different fuel ages and to align 

the flow vertically. These internal walls use a no-shear boundary condition since all flow 

resistance within the racks and fuel is accounted for by the porous model. 

Fresh air enters the building through the 4 m2 inlet vent (Figure 3). At this location, the flow is 

aligned normal to the building wall and has a turbulence level of 10%. The inlet velocity is set 

to yield a volume flow rate equivalent to 2 building volumes per hour (except as noted in 

sensitivity studies). The input air is 27 0C. The prescribed inlet velocity is 1.1 m/s. Air exits the 

model at the opposite end of the building. The exit vent is modeled as a pressure boundary set 

to atmospheric conditions (101325 Pa, 27 0C). 

4.3 Material Properties 

The working fluid is air modeled as an ideal gas. The temperature-dependent thermal 

properties are given in Table 2. These properties are obtained from White's Viscous Fluid 

FloW. The FLUENT code uses linear interpolation between points in the material property 

tables. For temperatures above the maximum tabulated value (1000 0C), the FLUENT code 

uses the thermal properties associated with the maximum value in the table. This assumption 

is valid as long as the temperatures do not go significantly above the values in the table. 

Table 2. Material Properties, Air (White, 1988) 

Temperature 'C Cp (J/kg-K) k (W/m-K) viscosity (Pa-s) 

0 1003 .0242 1.72 e-05 

100 1010 .031 2.17 e-05 

200 1024 .0376 2.57 e-05 

300 1045 .0439 2.93 e-05 

400 1069 .0497 3.25 e-05 

500 1093 .0553 3.55 e-05 

600 1114 .0609 3.83 e-05 

700 1135 .0659 4.09 e-05 

800 1153 .0703 4.34 e-05 

900 1170 .0747 4.57 e-05 

1000 1184 .0786 4.79 e-05

7
4.4 Fuel and Racks

The fuel and rack structure are not modeled in detail; the geometry of this region is too 

complex. CFD modeling of the details of every bundle in a fuel pool is impractical because of 

memory and processor limitations. To simplify the model, a homogenous nonisotropic porous 

medium is used to represent the racks and fuel. The volume occupied by the racks is replaced 

by the porous medium designed to affect the flowfield as the rack structure and fuel would. A 

porous material does not affect the flow area in the FLUENT model. Flow resistance is set 

through user-defined viscous and inertial flow loss coefficients. Details of the resistance 

coefficient determination are given in Appendix A. Lateral flow resistance in the porous medium 

is set two orders of magnitude higher than the vertical resistance to ensure vertical flow through 

the racks. Several vertical internal walls are also included in the racks to further guarantee a 

vertical flow path. The fuel is centered vertically in the racks and has a uniform active heated 

length of 3.8 m. The width of the active fuel region is 9.24 m, as shown in Figure 1. The 

individual fuel widths are defined to keep the fuel proportions the same as in Table 1. Fuel 

regions fl to f4 are 1.76, 0.5874, 0.5874, and 6.3052 m wide, respectively. Heat is added to 

the flow as a uniform volumetric heat generation term throughout the volume occupied by the 

fuel. Axial heating distributions are not considered. Radiation and clad chemistry are not 

included in the model. 

Thermal inertia and thermal conduction within the fuel and racks are accounted for in the 

FLUENT model through the porous media options. The material properties of the medium and 

the porosity of this region are prescribed to account for the thermal inertia of the racks and fuel 

and to approximate axial thermal conduction. The FLUENT porous media model assumes 

isotropic conductivity so the single value affects axial and lateral conduction. A sensitivity study 

on the thermal conductivity is done to give further insights into the significance of this 

parameter. Details of the porous media assumptions used to account for thermal inertia and 

thermal conductivity in the rack region are given in Appendix B. 

4.5 Turbulence Modeling 

The k-epsilon turbulence model is used for this analysis with standard wall functions. This 

model is selected because of its robustness and ease of application. It is acknowledged that 

the k-epsilon turbulence model is not tailored to this type of problem. Experimental data are 

needed to refine the turbulence modeling beyond this standard approach. This subject is 

beyond the scope of this initial study. 

The flow up (or down) through the bundles is expected to be laminar because of the low 

Reynolds number in this region. For example, if the velocity of 400 'C air is .243 m/s along the 

tube bundles (0.1 m/s in the open porous region of the FLUENT model), the Reynolds number 

based upon hydraulic diameter would be close to 50, indicating laminar flow. To account for 

this laminar flow, the porous media region representing the racks and fuel is modeled as a 

laminar region. In this region, turbulence production and turbulent viscosity are disabled. The 

FLUENT code model transports the turbulent properties through the region unchanged.

8
4.6 Solution Convergence

Final results from this analysis are somewhat unsteady. Variations of +/- 20 0C in the hottest 

fuel temperatures are typical. To obtain final results, the model is executed until the oscillations 

are in a steady-state pattern. At this point, variables are recorded during one or more 

oscillations and an average value is recorded. Initial convergence of the solution is monitored 

by observing temperatures and velocities at various points in the domain as a function of time. 

In addition, the residuals (a measure of error) are monitored and minimized for each solution. 

As a final check, the overall mass and energy balances are verified. 

4.7 Grid Independence 

A grid independence study is not feasible due to the physical scale of the model and the limited 

computer resources. A high quality mesh is used for this analysis to minimize any grid affects 

on the results. The effect of the assumptions on the final solution are expected to far outweigh 

any small effect of the grid on the predictions. Qualitative information on grid independence is 

given below for completeness. 

The final grid for this study contains over 675,000 finite volume cells. In the process of 

constructing this grid, over 10 separate grids are created and tested. The final grid is a 

compilation of the lessons learned during this process. The lessons learned included node 

density requirements for wall functions and large gradients. The final mesh is created with as 

many cubical elements as possible. Grid stretching and skew are minimized. Mesh size is 

reduced at the walls to accommodate the wall functions and transitions away from the wall are 

limited to growth rates between 5% and 20%. To transition between nearly cubical elements of 

different scales, tetrahedral regions are constructed with growth rates between 5% and 20%. 

The largest gradients are found just above the fuel racks. In this region, grid adaption is used 

twice to reduce the grid size to approximately 0.01 m. Unphysical results were observed in 

some of the earlier models without this adaption. Grid adaption is also used to refine the mesh 

in the open region around and below the racks. 

5 COMPLETED PREDICTIONS 

CFD results are most valuable when a series of predictions are made with varying assumptions 

and inputs to assess their impact on the predicted results. This is especially true when there is 

a lack of data available for model validation. Sensitivities of the final predictions to changes in 

assumptions and parameters gives valuable insights into the predictions. Table 3 gives the 

major parameters and assumptions varied for each prediction. As noted in Table 3, several 

sensitivity studies are completed which include variations of the following parameters: 

decay time 

ventilation rate 

flow resistance within the fuel and racks 

heat transfer coefficient for external building walls 

fuel loading pattern (location of hottest fuel) 

porous media thermal conductivity 

fuel burnup

9
Table 3. Completed Predictions

Decay Fuel Vent Overall h Sensitivity study notes 

time burnup (bldg./ (W/m 2-K) 

(yrs) (GWd/ hr) 

MTU) 

f2v2hO 2 40 2 0 2 yr, adiabatic walls 

f3v2hO 3 40 2 0 3 yr, adiabatic walls 

f4v2hO 4 40 2 0 4 yr, adiabatic walls (base case) 

f6v2hO 6 40 2 0 6 yr, adiabatic walls 

f4vl hO 4 40 1 0 vent 1 building per hour 

f4vlp5hO 4 40 1.5 0 vent 11/2 buildings per hour 

f4v2p5hO 4 40 2.5 0 vent 21⁄22 building per hour 

f4v2h2 4 40 2 2 wall heat transfer (h = 2 W/m2-K) 

f4v2h4 4 40 2 4 wall heat transfer (h = 4 W/m2-K) 

f4v2hOrp 4 40 2 0 rack/fuel flow resistance +20% 

f4v2h0rm 4 40 2 0 rack/fuel flow resistance -20% 

f4mv2h0 4 40 2 0 hottest fuel in center of pool 

f4v2hOkm2O 4 40 2 0 porous conductivity -20% 

f4v2hOkm5O 4 40 2 0 porous conductivity -50% 

f4v2hOb5O 4 50 2 0 fuel burnup = 50 GWd/MTU 

f2v2h2 2 40 2 2 2 yr best estimate, (h = 2 W/m 2-K) 

f3v2h2 3 40 2 2 3 yr best estimate, (h = 2 W/m 2-K) 

f6v2h2 6 40 2 2 6 yr best estimate, (h = 2 W/m2-K) 

The 4 year case, f4v2hO, is used as the base case. Most sensitivity studies begin with or pass 

through f4v2hO. Best estimate results use a wall heat transfer coefficient of 2 W/m 2

-K.

10
6 SUMMARY OF ASSUMPTIONS AND LIMITATIONS

The completed CFD predictions give valuable insights into the three-dimensional natural 

circulation air flows. These flows provide the primary heat removal mechanism in a postulated 

complete loss of spent fuel pool coolant accident. The final results, however, are still 

dependent upon the modeling assumptions and limitations. To ensure that these issues are 

considered along with the final results, some major assumptions and limitations are listed in 

Table 4. 

Table 4. Assumptions and Limitations 

Assumptions Notes/Limitations 

building and pool dimensions building and pool are within typical range but do 

not represent all configurations 

4200 fuel assemblies, high density racks no rack-to-rack spacing, no open areas in pool 

flow resistance in rack/fuel region determined from a typical configuration, does not 

apply to all rack designs 

0.15 m (6 in.) gap from rack to pool wall typical for a full pool; some gaps reported as 

small as 0.05 m (2 in.) 

no structures within building plant-specific structures could affect the flow 

model pool and containment above pool containment only modeled from pool surface up, 

surface level plant-specific features could affect solution 

2 buildings per hour ventilation rate typical of an operable ventilation system, location 

of vents could affect the flow patterns 

fuel arranged in contiguous regions by no checker boarding of fuel is attempted, no open 

age cells are in racks 

decay heat levels average value used for each fuel grouping, best 

estimate for expected high burnup fuel 

porous media rack structure is homogenized, heat is added 

directly to the flow, fuel surface is not modeled, 

flow resistance is modeled, isotropic conduction 

external pool walls are adiabatic only heat transfer from upper building and ceiling 

is considered 

uniform heat generation in fuel axial fuel peaking is ignored 

radiation neglected, can be important at elevated T 

clad oxidation reaction neglected, can be important at elevated T 

steady-state solution sought no transient information available during heatup 

CFD modeling assumptions and grid no data used for assessing turbulence modeling, 

grid independence is not verified

11
7 RESULTS

The base case, f4v2hO, is described in detail to illustrate the flow phenomena predicted and to 

define some of the parameters used to describe the results. Next, the results for each of the 

predictions are summarized in tables. Finally, the sensitivity study results are summarized. 

7.1 Base Case Results (f4v2hO) 

The base case uses a 4 year decay time with the hottest fuel located along the left edge of the 

pool. All external walls are adiabatic. The ventilation rate is based upon two building volumes 

per hour, which corresponds to an inlet air velocity of 1.1 m/s. The inlet air temperature is 27 

0C and the flow has a 10% turbulence level at the inlet. The inlet vent flow area is 4 M2 . 

Important parameters and phenomena predicted for this case are described below. 

7.1.1 Maximum Temperature 

The maximum temperature predicted during this scenario varies with time. The peak 

temperature predicted is 540 0C with variations of +/- 20 °C. Peak temperature is observed 

near the top of the hottest fuel and on the symmetry plane. Figure 5 shows the computational 

domain with temperature contours displayed on the surface of the active fuel region. The 

hottest regions are lightest (whitest) in the figure. 

540 

498 

456 

414 

372 

330 

1288 

246 

204 

162 LX 

120 

Figure 5. Temperature Contours on Active Fuel Region 

Walls between the fuel types are left in place as a visual aid. There is an additional vertical wall 

in fuel region fl and several additional walls in fuel region f4. To clarify the image, these walls 

are omitted. The hottest fuel, fl, is on the left, with a temperature approaching 540 0C. The air 

entering the bottom of the hottest fuel is approximately 125 °C. Actual fuel surface 

temperatures can be approximated using Newton's law of cooling. The Nusselt (Nu) number is 

approximately8 9. The difference between the predicted air temperature (assumed to be the 

mixed mean) and the surface temperature is predicted to be less than 10 0C. This difference is 

considered negligible and the maximum temperatures are therefore simply reported as the 

maximum predicted air temperatures.

12
7.1.2 Global Flow Pattern

The overall flow, although unsteady, displays a predictable pattern. Figures 6, 7, and 8 

illustrate the principal convective flow paths. Figure 6 shows the volume of air that is predicted 

to have temperatures in the range of 132 to 500 0C. These limits are chosen to highlight the 

hot plume of air rising from the side of the pool containing the hottest fuel. The rising plume 

stays near the left pool wall until it enters the main building. At this point, the hot air flows 

towards the roof of the containment with a slight preference for the outlet vent. The rising hot 

air hits the ceiling and spreads to cover the entire ceiling with a layer of hot air. Although the air 

temperature peaks at over 500 0C in the fuel racks, the rising plume of hot air cools rapidly as it 

rises. The maximum temperature is 200 °C just a few meters above the fuel racks. As the hot 

plume reaches the top of the pool, the maximum temperature is 160 0C. The maximum plume 

temperature is 135 0C when the plume hits the ceiling. Temperatures in the containment are 

predicted to stratify significantly. The temperature of the hot air on a horizontal plane 0.5 m 

below the ceiling ranges from 130 to 135 0C. Eliminating the small region occupied by the 

plume reduces the temperature variation on this horizontal plane to less than 2 °C. 

4.85e+02 4.569+02 

4.26e+02 

3.97e+02 

3.68e+02 

3.388+02 

3.09e+02 

2.79&+02 

2.509+02 

2.20e+02 Contours of Static Temperature (c) 

1.91 e+02 

1.61 e+02 

1.32e+02 

Figure 6. Temperature Contours Showing Rising Hot Plume 

Figure 7 shows the inlet air entering the containment and falling to the floor of the building. This 

relatively cool air spreads around the rising hot plume exiting the pool and falls into the pool on 

the opposite side. The air on the floor of the containment is relatively well mixed. The 

temperature range observed on a horizontal plane 0.5 m above the containment floor is 80 to 

160 0C. Eliminating the small regions where the inlet plume first hits the floor and where the hot 

plume rises through this plane results in a temperature range at this elevation of 90 to 110 0C. 

The average temperature difference between the floor and the ceiling in the containment is 

approximately 30 °C. 

A detailed view of the flowfield just above the containment floor is given in Figure 8. This figure 

shows the velocity vectors on a plane 0.5 m above the floor. The cold inlet plume is predicted 

to hit the floor and spread out in all directions. The flow travels across the containment floor

13
and around the rising plume. The flow falling into the pool is also clearly visible. The pool 

edges are not visible. The increase in vector density in the rectangular region around the pool 

is the result of a mesh density variation near the pool. The mesh density variation and the pool 

edges are visible in Figure 6. The cold air that enters the fuel pool hits the top of the racks (like 

a stagnation point) and then flows outward towards the pool walls. Most of the cooling air 

travels down around the racks (along walls) to the floor of the fuel pool and then up through the 

fuel. 

1.11e+02 

1.04e+02 

9.70e+01- - - -

8.98e+O1 

8.26e+-01 :T 

7.54 ÷ l....... . 

5.38e+01 

4.66e+o1 Contours of Static Temperature (c) 

3 F94g+01 

3.22e+01 

2,50e+01 

Figure 7. Temperature Contours Showing Cooling Air Flow

Figure 8. Velocity Vectors in a Plane 0.5 m Above the Containment Floor

14
7.1.3 Global Energy Balance

To verify a global energy balance, the first law of thermodynamics is applied to the system 

using a control volume approach. The outlet temperature is computed and compared to the 

prediction from the FLUENT code. The relevant energy terms are the energy input from the 

fuel, the energy lost at the external walls by convection, and the energy transported away by the 

ventilation system. In this case, all the heat generated by the fuel is carried away by the 

natural circulation flows since the walls are adiabatic. The first law of thermodynamics applied 

to this steady-state/steady-flow problem results in the following balance of energy equation: 

Q fuel - Q walls + m(hinti, - houtete) = 0 

The terms in the equation represent the energy from the fuel, the energy lost at the walls due to 

convection, and the net energy change between the inlet and the outlet flows. For this case, 

total heat from the fuel, Qfue, is 551,474 W. The wall heat transfer, Qe,,wa, is 0. Assuming 26.85 

'C air entering the building through a 4 m2 inlet vent at 1.10 m/s results in an inlet mass flow 

rate of 5.178 kg/s. From thermodynamic tables, the inlet enthalpy (h,,net) is 300.19 kJ/kg. Using 

these values for the total heat load, the mass flow rate, and the inlet enthalpy, the energy 

balance predicts an outlet enthalpy (ho,,tet) of 406.69 kJ/kg. From thermodynamic tables, this 

translates into an outlet temperature of 132.5 0C. The predicted outlet temperature from 

FLUENT is 132.7 +/- 5 'C. The FLUENT code maintains a global energy balance. 

7.1.4 Flow In and Around Racks and Fuel 

The flow from the upper building enters the pool and falls to the racks as shown in Figure 7. 

The flow turns at the rack surface. Some of the flow is entrained by the rising hot air and exits 

the pool. Some travels down and around the fuel racks to the pool floor, where it spreads out to 

cover the region below the racks. As air is pulled up through the fuel and racks, it forms the hot 

plume visible in Figure 6. The regions of up and down flow in the rack region are illustrated in 

Figure 9. This figure shows contours of vertical velocity on a horizontal plane at the top of the 

active fuel region (z = 4.2 m). 

0.12 upfl ow > O.12m/s downflow<-0.12t 

o.09 downcorner flows go beyond contour range 

o0o. 

0.04 

0.00 

-0.03 

-0.05 t 

-0.06 w- . 

-0.09 

-0.12 

Figure 9. Vertical Velocity (m/s), Horizontal Plane, Top of Active Fuel

15
The contour ranges are limited to +/- 0.12 m/s to highlight variations in the fuel region. 

Downcomer velocity magnitudes go beyond these limits, however. Regions where the velocity 

is above or below the contour ranges are displayed with the maximum or minimum contour 

range. The velocity in the fuel is positive (upward) everywhere except a small region near the 

right downcomer. 

A slight upflow is predicted for the left downcomer. This upward flow is the result of several 

factors. First, the strong plume formed by the hottest fuel blocks downflow on the left side of 

the pool. In addition, the left downcomer is subject to some buoyant forces due to heat transfer 

from the hottest fuel adjacent to this region. Entrainment of air in the rising plume on the left 

side of the pool also promotes upflow in the left downcomer. Finally, a strong right-to-left flow 

under the fuel racks feeds the flow up through this downcomer. 

Most of the flow to cool the fuel enters the region below the racks from the downcomer on the 

right side of the pool. This region is supplied by a significant downflow of cooling air from the 

upper building. The air entering the region below the racks spreads across the floor, as shown 

in Figure 10. In this figure velocity vectors are shaded with temperature on a plane midway 

between the floor of the pool and the lower edge of the racks (z = 0.075 m). Figure 10 clearly 

shows that the source of cooling air is the right downcomer. 

158,40 

. .. . . ..... ... . 

. 148.80 . .. ..... " .. e-. ..- 

145. 0 .. . ... .. ... ... "- .... - .- •... - ....... -.-. - + ...... -,..... <e .--.-. • -.. 4: ..... .. • 

. . ._-. - 4. 

9139.20 ... . . -- .. .-- E 136.00 .. . -. -• .. .. - -- .. 

-2 132.80 . -_ -- • . 4 .---. 4----- .... - (• .. . . . ÷ ~ 4..... . .... -- 4- 4..... . 43- •• - .•- - .--- •- t- .-.. 4-

(1) 1 2 9 .6 0 . - .- .- - , - - -- ..... .. 

126.-40- .4 

126...40 - * • - -....- 44*_ • • -___= . 4. - -4-- .r-- -4 , .....--- - , •--'--4- r3"" 

126.40 " -- 4.- 4 -- 4• ÷ ..... ....... * "---- -- " •

123120 ... S.. . - - --- -- -- - - - "-- '- " ..- . 4-- -4-- -- 4-- -- 4,-1 ,. -4---. 4-.-.-- --

120,00 

Figure 10. Velocity Vectors Shaded With Temperature Below the Fuel Racks 

Flow in the third downcomer region (top of Figure 9) is more complex. No strong up or down 

flow is observed. The left side has a slight upflow and the right side has a slight downflow. 

Large recirculations are observed in this region, making it difficult to generalize the flow path. 

All fuel and downcomer flows are summarized by the net mass flow rate through each region. 

Mass flow rates are reported on a horizontal plane at the top of the active fuel region (z = 4.2 

m), with positive values indicating upward flow.

16
Net mass flow through a horizontal plane at the top of the active fuel is given in Table 5 for the 

four fuel regions and three downcomer regions. The predicted mass flow rates are slightly 

unsteady. Reported values represent an average mass flow rate for a series of predicted 

values over time. 

Table 5. Mass Flow Rates at Top of Active Fuel, Base Case

7.1.5 Pressure Above and Below Rack 

Spatial pressure variations are predicted above and below the racks. These predictions can be 

used to assess the applicability of the constant pressure boundary condition used in other spent 

fuel pool models in these regions. All predictions are given as a gauge pressure referenced to 

the atmospheric pressure at the outlet vent. 

Figure 11 shows the predicted pressure on a horizontal plane at an elevation of z = 4.5 m. This 

plane is 0.05 m (2 in.) above the top of the racks. Pressures vary from 15.25 to 18.25 Pascal 

(Pa). The average pressure at this elevation is 16.7 Pa. The right side of Figure 11 shows the 

highest pressure near the symmetry plane. Figure 7 shows why: the relatively cool air flows 

down the right side of the pool and impinges on the racks there. This region of high pressure is 

the result of the stagnation region created by the falling plume as it turns 90 degrees to flow 

along the top of the racks. The stagnation point is on the axis of symmetry near the right 

downcomer. Flow vectors radiate from this point in all directions as the velocity vectors, not 

shown, clearly demonstrate. 

The magnitude of the pressure variation across the top of the racks is predicted to be 

approximately 3 Pa. This value is considered significant relative to the pressure differences 

encountered. For example, the predicted difference between the average

17

Region Mass Flow Rate Average Mass 

(+ upflow, kg/s) Flux (kg/m2 s) 

left downcomer 0.22 0.27 

right downcomer -1.58 -1.90 

rear downcomer -0.18 -0.13 

(top of Fig. 9) 

net downcomer >> -1.54 -0.51 

fl 0.41 0.043 

f2 0.13 0.041 

f3 0.12 0.037 

f4 0.89 0.018 

net fuel >> 1.54 0.031
18.12

17.86 

17.60 ! 

"17.34 

17.08 

16.82 

16.56 

16.30 

16.04 

15.78 

15.52 

15.26 

15.00 

Figure 11. Pressure Contours Above the Fuel Racks, Pa (z = 4.5 m) 

pressure above and below the racks for the base case is 5 Pa. Flow losses are also in this 

range. Consider air at 900 °C flowing through the racks and fuel at 0.1 m/s. The pressure drop 

(no gravity) is approximately 15 Pa based upon the coefficients in Appendix A. If the air is 400 

°C and flowing at 0.01 m/s, the expected pressure drop is close to 1 Pa. Considering these 

pressure differences, the 3 Pa variation across the top of the racks could be significant. 

Figure 12 shows the predicted pressure contours below the rack at an elevation of 0.1 m. This 

plane is 0.05 m from the bottom of the racks. Pressures vary from 7.25 to 15.25 Pa; the 

average pressure is 11.8 Pa. The right side of Figure 12 shows the highest pressure at the 

bottom of the right downcomer. This is the region where the strong downflow hits the floor and 

turns to flow along the bottom of the racks. The high velocities leaving this region result in low 

pressures to the left of the stagnation region. The static pressure slowly increases as the flow 

slows down on its way to the hottest fuel region at the left. As noted earlier, the reported 

pressure variations at this plane are significant with respect to other relevant pressure 

differences in this region. 

7.1.6 Temperatures Above and Below Rack 

The temperatures above and below the rack are given by Figures 13 and 14, respectively. 

These results are highlighted to provide a basis for assessing the constant temperature 

boundary condition assumed by some spent fuel codes. Temperatures are given on the same 

horizontal planes used to illustrate pressure in Figures 11 and 12.

18
19

14.93 

14.29 

13.65 

i13.01 

12.37 

11.73 

11.09 

10.45 

9.81 

9.17 

8,53 

7.89 

7.25 

Figure 12. Pressure Contours Below the Fuel Racks, Pa (z = 0.1 m) 

316 0C 

298 

280i 

244 

226 

-208 

190 

172 

154 

136 

118 

100 

Figure 13. Temperature Contours Above the Fuel Racks (z = 4.5 m)
177.6 Oc

172.8 

168.0 

163.2 

158.4 

153.6 

148,8 

144.0 

139.2 

134.4 

129.6 

S124.8 

120.0 

Figure 14. Temperature Contours Below the Fuel Racks (z = 0.1 m) 

Temperature contours above the fuel racks (Fig 13.) are consistent with the description of the 

flow paths given earlier. The stagnation region on the symmetry plane near the right side 

downcomer has the lowest temperature. This is a result of the cold plume falling into the pool 

region and stagnating on the racks at this point (see Figure 7). As the flow travels across the 

racks to the left, it mixes with hotter air exiting the racks. This increase in temperature is clearly 

seen in Figure 13. The rising hot plume on the left side of Figure 13 is the hottest region at this 

level. This area roughly covers the hottest fuel region (f 1). The average temperature above the 

racks at this level is 217 °C. This value is computed over a series of iterations to get a 

representative value. The spatial variation in temperature at this level is roughly 200 degrees. 

Figure 14 shows the temperature contours below the fuel racks. The average temperature for 

this level is 137 °C. The spatial variation in temperature is nearly 60 degrees. Most of this 

spatial variation is accounted for by the hot spot near the top right hand corner of Figure 14. 

This region is caused by a small downflow predicted to occur in the fuel racks in the cells just 

above this region. This downflow is heated by the fuel before it enters the region below the 

racks. 

The temperature below the racks is nearly equal to the containment outlet vent temperature. 

This feature of these predictions is coincidental. There is a 30 degree average temperature 

difference between the floor and ceiling in the containment. The outlet vent is near the ceiling 

and is exposed to the hotter layers of air in the containment. The air which cools the fuel 

comes mainly from the floor of the containment. As it travels to the region below the racks, it 

rises in temperature to nearly equal the outlet vent temperature.

20
7.2 Tabulated Results for Each Case

The predicted values for each of the case studies described in Table 3 are outlined in Tables 6, 

7, and 8. Table 6 lists relevant temperatures and pressures. Table 7 summarizes mass flow 

rates through the fuel region. Table 8 is a summary of the global energy balance for each 

case. Each of the tabulated results represents an average value. The solutions are predicted 

to be naturally unsteady. Averages are obtained from approximately 1000 iterations recorded 

once the solution reaches a steady state pattern. The magnitude of the variation in depends 

upon the variable considered. The results for each sensitivity study are discussed following the 

tables. 

Table 6. Predicted Temperatures and Pressures 

case ID max T T at Average T (°C) Average P (Pascal) 

(°C) outlet 

(°C) above rack below rack above rack below rack 

f2v2hO 1047 176 297 177 18.9 15.2 

f3v2hO 686" 148 245 149 17.5 13.3 

f4v2hO 537 133 217 137 16.7 11.8 

f6v2hO 445 120 193 128 16.1 10.7 

f4vl hO 954 237 309 259 6.5 2.9 

f4vlp5hO 642 167 254 171 11.1 7.3 

f4v2p5hO 518 112 186 123 23.1 16.6 

f4v2h2 527 102 174 113 14.9 9.5 

f4v2h4 477 91 171 111 14.0 6.7 

f4v2hOrp 586 133 218 135 16.8 12.1 

f4v2hOrm 519 133 216 139 16.7 11.4 

f4mv2h0 567 133 212 138 16.4 12.1 

f4v2hOkm2O 543 133 215 137 16.7 11.6 

f4v2hOkm5O 570 133 207 140 16.7 11.9 

f4v2hOb5O 712 160 268 163 18.0 13.5 

f2v2h2 847 134 243 143 16.5 13.0 

f3v2h2 595 114 202 125 15.5 11.5 

f6v2h2 392 94 167 100 14.4 9.1

21
Table 7. Predicted Mass Flow Rates

case ID inlet flow rate through fuels (kg/s) downcomer flow (kg/s) 

(kg/s) 

fi f2 f3 f4 net left rear right net 

f2v2hO 5.18 .35 .11 .11 .81 1.39 .22 -. 14 -1.46 -1.38 

f3v2hO 5.18 .40 .13 .12 .90 1.55 .25 -.26 -1.53 -1.54 

f4v2hO 5.18 .41 .13 .12 .89 1.54 .22 -. 18 -1.58 -1.54 

f6v2hO 5.18 .39 .13 .12 .92 1.56 .22 -.20 -1.57 -1.56 

f4vl hO 2.59 .20 .07 .07 .74 1.09 -.99 -.33 .23 -1.08 

f4vlp5hO 3.88 .36 .11 .11 .84 1.42 .17 -.17 -1.41 -1.41 

f4v2p5h0 6.47 .42 .13 .12 .87 1.54 .31 -.2 -1.65 -1.54 

f4v2h2 5.18 .44 .14 .13 1.09 1.79 .17 -.48 -1.48 -1.79 

f4v2h4 5.18 .45 .14 .13 .93 1.64 .33 -.31 -1.65 -1.63 

f4v2hOrp 5.18 .36 .12 .11 .83 1.42 .19 -.05 -1.55 -1.41 

f4v2hOrm 5.18 .45 .14 .13 .96 1.68 .31 -.40 -1.58 -1.68 

f4mv2h_ 5.18 j.38 j.13 I 10 1.04 1.65 I.04 -. 19 -1.50 -1.65 

f4v2hOkm2O 5.18 .40 .13 .12 .88 1.53 .26 -.23 -1.55 -1.52 

f4v2hOkm5O 5.18 .40 .13 .12 .90 1.55 .30 -.37 -1.47 -1.55 

f4v2h~b50 5.18 .37 j.12 .11 .89 1.45 J.18 -. 13 I -1.53 -1.48 

f2v2h2 5.18 .46 .15 .14 1.12 1.87 .30 -.87 -1.28 -1.86 

f3v2h2 5.18 .49 .16 .15 1.18 1.97 .34 -.97 -1.33 -1.97 

f6v2h2 5.18 .45 .14 .14 1.09 1.82 .19 -.28 -1.72 -1.81

22
Table 8. Predicted Global Energy Balance Terms

case ID Heat Load from Fuel (kW) Convected Heat transfer 

fl f2 f3 f4 net enthalpy at walls 

f2v2hO 368.0 71.6 50.9 288.0 778.4 -778.4 0 

f3v2hO 248.0 55.9 45.5 283.7 633.2 -633.2 0 

f4v2hO 184.8 48.1 40.6 278.0 551.5 -551.5 0 

f6v2h0 136.4 39.5 36.3 272.2 484.5 -484.5 0 

f4vl hO 184.8 48.1 40.6 278.0 551.5 -551.5 0 

f4v1 p5hO 184.8 48.1 40.6 278.0 551.5 -551.5 0 

f4v2p5hO 184.8 48.1 40.6 278.0 551.5 -551.5 0 

f4v2h2 184.8 48.1 40.6 278.0 551.5 -395.1 -156.4 

f4v2h4 184.8 48.1 40.6 278.0 551.5 -338.9 -212.6 

f4v2hOrp 184.8 48.1 40.6 278.0 551.5 -551.5 0 

f4v2h0rm 184.8 48.1 40.6 278.0 551.5 -551.5 0 

[_f4mv2hO 184.8 1_48.1__J40.6 278.0 551.5_ -551.5 0 

f4v2hOkm2O 184.8 48.1 40.6 278.0 551.5 -551.5 0 

f4v2hOkm5O 184.8 48.1 40.6 278.0 551.5 -551.5 0 

f4v2hOb5O 230.4 61.0 J 52.1 352.6 696.0 -696.0 0 

f2v2h2 368.0 71.6 50.9 288.0 778.4 -563.2 -215.2 

f3v2h2 248.0 55.9 45.5 283.7 633.2 -458.6 -174.6 

f6v2h2 136.4 39.5 36.3 272.2 484.5 -352.2 -132.3 

7.3 Decay Time Sensitivity, Adiabatic Walls 

The primary sensitivity study focuses on the decay time, defined as the time since reactor 

shutdown. As decay time increases, the heat load from the fuel decreases. Decay times of 2, 

3, 4, and 6 years are studied. The physical and numerical conditions are identical for each 

case. Only the decay time is varied. Each case is based upon an assumed burnup of 40 

GWd/MTU. Predicted temperatures are plotted in Figure 15. Spline fits connect the data. The 

maximum temperature decreases significantly (511 0C) between 2 and 4 years. From 4 to 6 

years, the drop in temperature is only 91 0C. The shape of this curve closely resembles the 

decay heat curve. It is important to keep in mind the limitations (noted earlier) of these 

predictions at elevated temperatures. The 2 and 3 year results reach temperatures where the 

lack of radiation and chemistry models is expected to affect the solution.

23
1100

-- maximum T 

900 -- outlet vent T 

90.- average T above racks 

-7-: average T below racks 

S700 

500 sp line fits 

( • outlet T and T below racks 

300 

. 

-U -,,• -- -

100 ...... 

2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 

decay time (years) 

Figure 15. Predicted Temperatures vs. Decay Time, Adiabatic Walls 

Temperatures above and below the rack differ by nearly 80 0C. These temperatures decline 

with decay time as expected. The rate of decline also decreases significantly after 4 years. 

The outlet vent temperature is predicted to decline in the same manner. One unexpected 

prediction is the near match between the outlet temperature and the average temperature 

below the racks. This coincidental result is discussed in a previous section. 

The average pressure above and below the racks shows a slight increase as the decay time 

decreases. This is mainly because the pressure drop at the outlet vent increases as the decay 

time decreases. The outlet vent is at the end of a 4 m convergent duct protruding out from the 

containment wall. The pressure drop in this duct increases with the increase in viscosity and 

velocity associated with lower decay times. The increase in flow resistance through the outlet 

duct slightly increases the global containment pressure. 

The flow patterns for the 2, 3, 4, and 6 year cases are qualitatively the same. Downcomer 

flows are predicted to distribute in a similar manner. The left downcomer has a slight upflow for 

each case. The rear downcomer has a similar net downflow. Cooling flow enters the region 

below the racks mainly through the right downcomer. No clear trend with decay time is 

observed. The prediction of less total mass flow for the 2 year case is unexpected. The 2 year 

case does show larger velocities in the fuel region. The decrease in density, however, results 

in less total mass flow. The increased temperature associated with the 2 year case results in 

an increase in viscosity. The increased viscosity and velocity in this case increases the 

pressure-loss terms for flow through the racks. This factor is assumed to be partly responsible 

for the decreased mass flow associated with the 2 year case.

24
Flow through each of the fuel regions shows a similar pattern. No clear trend is observed with 

decay time variation. The noticeable exception is the reduced flow rate predicted for the 2 year 

case. Predictions for the 3, 4, and 6 year cases show nearly identical flow rates through each 

of the fuel types. The similar flow rates indicate that the increases in velocity for the lower 

decay times (higher T) are offset by similar decreases in density. 

7.4 Ventilation Rate Sensitivity 

An important input parameter is the ventilation rate. The base ventilation rate assumes 2 

building volumes of cool air enter the building each hour. Additional predictions are made with 

1, 1.5, and 2.5 building volumes per hour. Each case assumes a 4 year decay time and all 

external walls are adiabatic. Predicted temperatures are plotted in Figure 16.

1100

M 

0 

E 

-D 

ci 

ci 

L._

900 

700

500 4

300 -

100

1.0 1.5 2.0 2.5

ventilation rate (building volumes per hour) 

Figure 16. Sensitivity of Temperatures to the Ventilation Rate 

The maximum temperature is significantly affected by the ventilation rate. The largest affect 

occurs when the ventilation rate is reduced from 1.5 to 1 building volume per hour. This trend 

suggests that reducing the ventilation rate below the 1 building per hour rate would significantly 

increase the temperature. This is not done because the model is not accurate at higher 

temperatures. The difference in the maximum predicted temperature between 2 and 2.5 

buildings per hour is 19 0C. This is nearly equal to the difference in the predicted exit vent 

temperature for these two cases. It appears that the effect of increasing the ventilation rate 

beyond 2 buildings per hour can be approximated by globally adjusting the temperatures to the 

change in the outlet temperature determined from a global energy balance.

25

4 year decay time, 40 GWd/MTU 

..- maximum T 

-- outlet vent T 

- average T above racks 

..-- average T below racks 

outlet T and 

T below racks
The pressure in the containment increases slightly as the flow rate increases due to increased 

pressure losses in the outlet vent. The difference in the average pressure between the upper 

and lower racks also increases with flow rate. This is attributed to the increased mass flow 

through the downcomer region and the associated increased pressure drop. 

At the lower ventilation rate, the mass flow through the fuel drops off. The increased 

temperature reduces the density and this affects the mass flow. As the ventilation is increased, 

the mass flow through the fuel increases up to a point. No difference is observed in the mass 

flow rate through the fuel with the 2 and 2.5 building per hour ventilation rates. It appears that 

ventilation rates beyond 2 building volumes per hour don't significantly change the mass flow 

rates through the fuel or the net fuel temperature rise. The net fuel temperature rise is defined 

as the difference between the outlet vent temperature and the peak fuel temperature. Outlet 

vent temperature is considered a measure of the overall containment temperature. 

7.5 Wall Heat Loss Sensitivity 

A sensitivity study is done on the overall heat transfer coefficient between the containment and 

the exterior environment. Containment wall and ceiling heat transfer coefficients of 2 and 4 

W/m2

-K are applied to the base 4 year case. Predicted temperature results are plotted in 

Figure 17. The temperatures drop predictably as the heat transfer coefficient is increased. The 

drop in the maximum predicted temperature is most significant when the coefficient is changed 

from 2 to 4 W/m 2

-K. The temperatures above and below the racks respond differently. These 

temperatures change most dramatically when the heat transfer coefficient is changed from 0 to 

2 W/m 2-K. The global energy balance also changes most when the coefficient is changed from 

0 to 2 W/m 2-K. 

600 

- 500 

4 year decaytime, 40 GWd/M TU 40o0. 

maximum T 

E outlet vent T 

-U-i average T above racks 

- average T below racks 

a_) 

L-~ -0 

100 

13• 00 . .. . -- .. .. . .. - . - . - . --.. .. .... .. -

0 1 2 3 4 

containment wall and ceiling heat transfer coefficient (W/m 2-K) 

Figure 17. Sensitivity of Temperatures to Overall Heat Transfer Coefficient

26
The mass flow rates through the fuel region increase slightly as the heat transfer coefficient is 

increased. The limited data show no clear trend. The increase in the mass flow rate through 

the fuel is attributable to the lower containment temperatures resulting from the wall heat 

losses. The lower temperatures produce stronger natural circulation flows, which aid the fuel 

cooling. 

The sensitivity of wall heat loss is also observed by comparing the adiabatic decay time 

sensitivity results (f2v2hO, f3v2hO, f4v2hO, f6v2hO) with the best estimate results (f2v2h2, 

f3v2h2, f4v2h2, f6v2h2). These models differ only by the wall and ceiling heat transfer 

coefficient which is set to 2 W/m 2-K in the best estimate predictions. The data are given in 

tables 6, 7, and 8. The most significant difference is the drop in peak temperature when the 

heat transfer coefficient is applied. This temperature reduction is greatest at the 2 year decay 

time due to the higher containment temperatures in this case. The peak temperature of the fuel 

drops by 200 0C when wall heat transfer is applied using fuel with a 2 year decay time. At 6 

years decay time, the peak temperature drops by 53 0C when wall heat transfer is turned on. 

Temperatures at the outlet and above the fuel racks respond in a similar manner. No clear 

trend is observed in the average temperature computed below the fuel racks. Mass flow rates 

through the fuel and downcomers increase in the cases with wall heat transfer. This increased 

mass flow and reduced temperatures in the containment and fuel region account for the 

reduced peak fuel temperatures predicted. 

7.6 Rack and Fuel Flow Resistance Sensitivity 

The racks and fuel are modeled as a porous medium. The complexity of the physical structure 

in this area makes a detailed model impractical. A porous resistance is put into this region of 

the model to account for the viscous and inertial flow losses. Porous resistance is considered a 

significant parameter in these predictions. Details of the determination of the resistance 

coefficients are given in Appendix A. There is some uncertainty in the data used to determine 

the loss coefficients, and there are a variety of rack and fuel bundle designs that affect the flow 

loss. This sensitivity study shows the effect of flow loss variations on the predicted 

temperatures. 

Using the base case as a starting point, two additional predictions are done with the resistance 

coefficients changed by +/- 20%. The increased resistance reduces the mass flow rate which 

increases the temperature rise of the air flowing through the racks. Similarly, with less 

resistance the mass flow is increased and the temperature rise is reduced as the data in Table 

7 demonstrate. Figure 18 shows the predicted temperature for these cases. The maximum 

temperature increases with the increased resistance, as expected. This trend is attributable to 

the mass flow rates discussed above. The largest change occurs when the resistance is 

increased by 20%. The building outlet temperature is unaffected since the overall energy 

balance is unchanged. The average temperatures above and below the racks show no 

significant change.

27
600

500 

4 year decay time, 40 GWd/M TU 

0 400 

C--- -u - -- -- 3- 00 - vrg- eo ak 

100 

80% 100% 120% 

resistance coefficients 

Figure 18. Sensitivity of Temperatures to Flow Resistance 

7.7 Hottest Fuel Location Sensitivity 

The hottest fuel is placed along the left wall of the fuel pool for each of the predictions 

described above. The hottest fuel represents the last fuel taken from the reactor when the 

reactor is shut down. The fuel pool is assumed to be filled from right to left (Figure 1), so that 

the last fuel taken from the reactor is on the left. As a sensitivity study, the hottest fuel is 

moved to the center of the pool. Fuel region f2 is next to the hottest fuel on the left. Fuel 

region f3 is to the right of the hottest fuel. The remaining rack locations are filled in with the 

oldest fuel, f4. The total fuel load and quantities of each type of fuel are unchanged by this 

arrangement. 

The most conservative (hardest to cool) location for the hottest fuel is not obvious. Placing the 

hottest fuel in the center of the pool puts the fuel farthest from the downcomer regions 

supplying the cooling air flow, which could make it harder to cool. However, placing the hottest 

fuel next to the downcomer can block the cooling air from flowing down through this 

downcomer. This blockage can result from the rising hot plume. There is no benefit to putting 

the hot fuel next to a downcomer in this case. Two fuel positions are compared here to show 

the effect of fuel position on the predicted results. 

The base case assumes the hot fuel is along the left wall, as described above. Downflow to the 

left downcomer is blocked off by the rising plume. A small upward flow is established in this 

downcomer by entrainment and other factors. Cooling air from the right downcomer travels 

along the pool floor under the racks to reach the hottest fuel (see Figure 10). The average 

peak temperature in the fuel region is 537 'C.

28
The predictions done with the hottest fuel in the center of the pool are very similar to the base 

case predictions described above. The hot plume that rises from the center of the pool is swept 

to the left wall of the pool before rising into containment. The pattern is similar to Figure 6. 

Cold inlet air falls to the floor of the containment and travels around the pool, as shown in 

Figure 8. The average maximum temperature in the fuel region for this case is 567 °C. The 

mass flow rate through the hottest fuel is slightly smaller when this fuel is in the center of the 

pool. The inlet flow path apparently affects the flow patterns in the fuel pool. 

7.8 Rack and Fuel Thermal Conductivity Sensitivity 

The racks and fuel are modeled as a porous medium because of the physical structure in this 

area is too complex to model. This porous model requires the user to specify properties to 

account for given thermal behavior. The thermal property assumptions are given in Appendix 

B. The thermal properties of the porous medium affect the thermal inertia and the conduction 

flux in the porous region. Thermal inertia is only a factor in transient terms so it is not expected 

to affect the steady-state results. The thermal conduction term, however, does move fuel 

energy away from the hottest regions. Thermal conductivity is approximated by assuming one

dimensional conduction in the vertical direction. The FLUENT code assumes isotropic thermal 

conductivity. The isotropic assumption used by the code and the one-dimensional assumption 

used to approximate the thermal conductivity are both questionable assumptions. The base 

case thermal conductivity is reduced by 20% and 50% in separate predictions to quantify the 

sensitivity of the results to this input parameter. 

The maximum temperature increases slightly as the thermal conductivity of the porous medium 

decreases (See Table 6). The reduction in thermal conductivity slows the conduction of heat 

away from the hottest regions and accounts for this increase. The maximum predicted 

temperature increases by less than 30 °C as the thermal conductivity in this region is reduced 

by 50%. This is considered a weak sensitivity. The global flow patterns, the integrated mass 

flow through the fuel regions, and the global energy balance are not significantly affected by 

changing the porous thermal conductivity. 

7.9 Fuel Burnup Sensitivity 

All predictions are done assuming an average fuel burnup of 40 GWd/MTU. A sensitivity study 

is done using an average fuel burnup of 50 GWd/MTU and a decay time of 4 years. The higher 

fuel burnup assumption, when applied in this model, increases the energy generation term in 

the fuel region. The energy associated with a fuel burnup of 50 GWd/MTU is summarized in 

Table 1. The total pool energy at 4 years for the 50 GWd/MTU burnup falls between the 2 and 

3 year cases with a burnup of 40 GWd/MTU. 

The predicted temperatures and pressures reported in Table 6 for the higher burnup case fall 

between the 2 and 3 year cases, which use a burnup of 40 GWd/MTU. This is consistent with 

the observation noted above concerning the total pool energy. The mass flow rate predictions 

from Table 7 also generally fall between the 2 and 3 year cases. Similar to the other cases, the 

majority of the cooling flow travels down the right hand downcomer to the region below the 

racks. The average maximum temperature in the fuel region for this case is 712 OC.

29
7.10 Best-Estimate Critical Decay Time 

Predictions are made at each of the four decay times using best-estimate parameters. These 

parameters include an overall wall heat transfer coefficient of 2 W/m 2-K and a ventilation rate of 

2 building volumes per hour. These results represent the best estimate of the fuel 

temperatures and the convective flowfields, given the assumptions and limitations of this model. 

The difference between these results and the decay time sensitivity results (outlined in section 

7.3) is the heat transfer coefficient on the building walls and ceiling. The heat loss from the 

walls and ceiling lowers the temperatures compared to the adiabatic decay time sensitivity 

outlined in section 7.3 of this report. The general features of the solution remain unchanged. 

The primary result is the decay time sensitivity given in Figure 19. The data from cases f2v2h2, 

f3v2h2, f4v2h2, and f6v2h2 are plotted. Critical decay times of 26 and 35 months are noted on 

the figure for temperature limits of 800 0C and 600 °C, respectively. Critical decay time is 

defined as the post shutdown time required to ensure fuel temperatures don't rise above 

predefined temperature limits after a complete loss of fuel pool coolant. The 26 month decay 

time, associated with the 800 0C temperature limit, has significant uncertainty due to the model 

limitations at elevated temperatures discussed earlier. More details of these results are 

discussed earlier in this report and the principal data are outlined in Tables 6, 7, and 8. 

1100 

.--- maximum T 

900 - outlet vent T 

0 o, goo gees ---- average T above racks D- - average T below racks 

500 ., 

C2. 300 

100

2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 

decay time (years) 

Figure 19. Best-Estimate Results for Maximum Temperature vs. Decay Time 

8 SUMMARY 

In support of the NRC rulemaking activity related to decommissioning, the Office of Research is 

using the FLUENT CFD code to study the air cooling phenomena of a spent fuel storage pool 

after a complete pool drainage. This study provides a solution for the critical decay time of a 

BWR spent fuel pool configuration. A series of sensitivity studies are completed to show the 

significance of some of the assumptions and parameters used to complete the predictions.

30
This study focuses on the potential modeling limitations of previous spent fuel heatup 

predictions which relied on simplified flowfield models. This three-dimensional CFD treatment 

avoids the need for many flowfield assumptions in the fuel and rack region. These results give 

significant insights into the spent fuel pool cooling issue. The natural circulation flows predicted 

are the most significant factor in removing heat from the fuel under these conditions. These 

results help reduce the uncertainty associated with this most significant aspect (natural 

circulation flows) of spent fuel pool cooling. However, the model does not include terms for the 

radiation and exothermic cladding reactions. These terms are significant at temperatures 

greater than approximately 600 0C. The present predictions form one piece of an overall 

solution to this issue. When using these predictions, the assumptions and limitations must be 

considered. 

The series of predictions completed form a database of results which help define the natural 

circulation flows associated with this spent fuel cooling scenario. These results provide a basis 

for determining the applicability of the simplified flowfield assumptions used in codes such as 

COBRA-SFS, SFUEL, and SHARP. The simplified constant temperature and pressure 

boundary conditions used in previous studies can be assessed using the quantitative pressure 

and temperature variations predicted in this study. The basis for the boundary temperatures 

and pressures applied in these codes can also be re-considered in light of the results obtained 

here. Sensitivity studies using these codes can be guided by the results and sensitivities of this 

current study. Sensitivity studies and analysis with these simplified flowfield models are beyond 

the scope of this report which focuses solely on the CFD model and predictions. 

A similar global flow pattern is predicted for most of the solutions. The primary flow travels 

down around the fuel racks and then up through the hot fuel. A hot plume rises from the pool to 

the containment ceiling and exits through the outlet vent. Temperatures, pressures, velocities, 

and mass flow rates are recorded for a series of cases designed to quantify the effect of 

important input parameters. Quantitative results are tabulated and compared. The peak 

temperature is the predicted result of greatest interest. 

The primary sensitivity study focused on the critical decay time. The maximum predicted 

temperature drops with decay time along a curve similar to a decay curve. Using adiabatic 

containment walls, the maximum predicted temperature is 800 0C at approximately 32 months 

decay time. Maximum predicted temperature is 600 0C at 40 months. Using the best-estimate 

wall heat transfer coefficient on the containment walls and ceiling, the critical decay times for 

the 800 and 600 0C temperature limits are 26 and 35 months, respectively. The 800 and 600 0C 

limits are reference points only. Determination of appropriate critical temperature limits is 

beyond the scope of this report. 

Containment ventilation significantly affects the predictions. Increasing ventilation rates beyond 

the baseline of 2 building volumes per hour has an insignificant effect. Decreasing the 

ventilation rate to 1 building volume per hour increases the predicted peak temperature by over 

400 0C. For a given fuel decay and burnup, the ventilation rate is the most significant 

parameter in this study. These results are completed with adiabatic containment walls. 

The heat transfer coefficient on the containment wall is varied from 0 to 4 W/m2-K. A value of 2 

W/m 2-K is used for best-estimate predictions. The variations in the heat transfer coefficient at 

48 months decay time causes the peak predicted temperature to drop by 60 0C as the overall

31
heat transfer coefficient is increased from 0 to 4 W/m 2-K. A larger temperature change is 

predicted at shorter decay times, where the containment temperatures are higher. 

Porous resistance and thermal conductivity are parameters computed by hand for use in the 

model. These parameters are varied to quantify the effect of their uncertainties on the final 

predictions. Varying the flow resistance in the racks up and down by 20% resulted in 

temperature changes of 49 and -18 0C, respectively. Reducing the expected thermal 

conduction parameter in the fuel racks by 50% increased the peak predicted temperature by 33 

'C. These variations are minor compared to the effect of other input parameters in the model. 

The decay curve used in this study assumed an average fuel burnup of 40 GWd/MTU. A 

sensitivity study completed with 50 GWd/MTU at 4 years decay time increased the peak 

temperature by 175 0C. The increased burnup significantly increases fuel power. The increase 

in temperature associated with the higher burnup fuel is consistent with the increased fuel heat 

load. Fuel burnup is a significant parameter in these predictions. 

9 REFERENCES 

1. Nourbakhsh, H.P., et al, "Analysis of Spent Fuel Heatup Following Loss of Water in a 

Spent Fuel Pool, A Users Manual for the Computer Code SHARP," NUREG/CR-6441, 

May 1998. 

2. Benjamin, A.S., et al., "Spent Fuel Heatup Following Loss of Water During Storage," 

NUREG/CR-0649, March 1979. 

3. Sailor, V.L., et al., "Severe Accidents in Spent Fuel Pools in Support of Generic Safety 

Issue 82," NUREG/CR-4982, July 1987. 

4. Fluent Inc., "Fluent 5 Users Guide," Canterra Resource Park, 10 Cavendish Ct, 

Lebanon, NH 03766, 1998. 

5. Anderson, J.D., Computational Fluid Dynamics: The Basics with Applications, 

McGraw-Hill Series in Mechanical Engineering, McGraw Hill, 1995. 

6. ASHRAE, "Handbook of Fundamentals', 1981. 

7. White, F.M., Viscous Fluid Flow, Mc Graw Hill, 1991. 

8. Todreas, N.E., Kazimi, M.S., Nuclear Systems I: Thermal Hydraulic Fundamentals, 

Hemisphere Publishing Corporation, 1990.

32
APPENDIX A 

LOSS COEFFICIENTS FOR FLOW THROUGH 

BWR FUEL STORAGE RACKS 

A.1 BACKGROUND 

The computational fluid dynamics (CFD) analysis of spent fuel pool cooling documented in this 

report uses a porous medium to represent the racks and fuel. The porous model allows the 

user to input flow resistance coefficients to represent viscous and inertial flow losses in the 

region. The porous medium represents the effect of the racks and fuel on the flow, without the 

need to represent the precise geometry of the fuel and racks. This analysis is done with the 

FLUENT CFD. FLUENT's porous model adds a source term to the momentum equations to 

account for viscous and inertial flow losses in the porous regions. The physical flow area is not 

changed so the fluid does not accelerate. Equation A-1 1 gives the FLUENT source term in one 

dimension. The viscous term is first, followed by the inertial loss term. 

dP dx = D .,uV + C .Y2 PV I VI eq A-1 

In three dimensions, the coefficients D and C represent matrices. These values are input 

parameters for the FLUENT porous media model. This appendix gives background information 

and assumptions used in determining D and C in the form required for the FLUENT model. 

A.2 GEOMETRY 

Fuel bundles are stored in the spent fuel pool in racks that > > 

keep the bundles upright and provide neutron absorption. ,I I 

In this analysis, a high-density rack is filled with 9x9 BWR > > ,1 >1 

spent fuel bundles. The bundle characteristics are typical 

for a BWR. 

Figure A-1 shows a BWR high-density rack. The square 

boxes composing the rack are 0.154 m (6.06 in.) from 

center to center. Rack height (w/o legs) is assumed to be 

4.3 m (14.1 ft). The active fuel length is 3.8 m (12.45 ft) 

and is assumed to be centered in the rack from top to 

bottom. Rack distance from the floor is 0.15 m (5.9 in.). 

The walls separating the fuel cells are solid. For the 

purposes of the CFD model, multiple racks are joined 

together to form a single rack of the type illustrated in 

Figure A-1. This single rack spans the entire pool. A gap of 

0.15 m (5.9 in.) is modeled between the rack and each side 

wall of the pool. Every rack cell contains a complete 9x9 

fuel bundle. The racks form a homogeneous nonisotropic 

structure whose properties can be defined by considering a 

single rack cell. 

Figure A-i. High-Density Rack

A-1
Figure A-2 shows the cross 

section of a 9x9 fuel assembly 

in a fuel rack. The dark regions 

on each side are boral inserts in 

the rack structure. The rack 

cell pitch is 0.154 m (6.06 in.). 

The fuel bundle is contained in 

a channel box with an inner 

dimension of 0.134 m (5.3 in.). 

There is a gap between the 

rack walls and the fuel channel 

box. For the purposes of this 

analysis, all fuel rods are 

considered to be full length. 

The rod OD is 0.1118 m (.44 

in.). Rod pitch is 0.1438 m 

(.566 in.). There are a total of 

74 fuel rods and 2 water rods in 

the 9x9 cross-section. Table A

1 gives other values used to 

describe the geometry for the 

purposes of determining the 

loss coefficients. The left 

column in Table A-1 gives a 

convenient Identification for 

each value.

>1 ___________ 3{.cy.z

AOt 

c20 OD-Ti-sq.£c' 0 ,rm 

(.9s ,,)

NOT TO SCALE 

I I 

Ij'ý L*- M CELL ID 

Id.- (.'5,9;9 SN) • 1' .9?. wAZ L mncE SPACING ,

Figure A-2. Cross-Section of Rack and Bundle

Table A-1. Values Used to Determine Loss Coefficient in One Rack Cell

D1 Flow area available to FLUENT model (0.153924)2 = 0.023692597 m2 

D2 Channel box flow area (no rods) (0.1340612)2 = 0.017972405 m2 

D3 Channel box perimeter 4 x 0.1340612 = .5362448 m 

D4 Cross-section area of rods (74 * rr d2

/4) 0.007259294 m2 

D5 Cross-section area of water rods 0.000973283 m2 

D6 Perimeter of all rods (74 * Tr d) + (2 * TT dj) 2.75457357 m 

D7 Total wetted perimeter (D3 + D6) 3.2908 m 

D8 Total flow area (D2 - D4 - D5) 0.0097398 m2 

D9 Hydraulic diameter (4 A / P)= 4*D8/D7 0.0118388 m 

D10 Flow area ratio at grid spacers 0.6 (area at spacers / area at tubes) 

D11 Flow area ratio at tie plates 0.5 (area at tie plate / area at tubes) 

D12 Lower orifice diameter (bottom of bundle) 0.09144m (3.6 in.)

A-2
A.3 LOSS COEFFICIENT DETERMINATION

The FLUENT code uses the loss coefficients, C and D, in equation A-1 to determine the 

momentum sink term used in the model. The velocity in equation A-1 represents the FLUENT 

solver velocity, which is the flow through the unobstructed flow area (D1 from table A-i). In 

reality, the velocity is higher for a given mass flow rate because of the reduced flow area of the 

physical structure (D8 in Table A-i). Figure A-3 shows a fuel bundle in a rack cell from the 

side. Viscous loss is accounted for along the tubes and side wall. The inertial losses are 

accounted for at the lower base hole, the lower and upper tie plates, and the seven grid 

spacers. Each of these losses is discussed below. 

A.3.1 Viscous Loss 

The viscous losses come primarily from flow along 

the tubes. It is assumed that viscous losses are 

uniform for the entire length of the racks. These'

losses are based on correlations for tube flow. . , 

Equation A-2 defines the pressure loss over a 

distance L along tube banks with hydraulic diameter d. -15413S sp/er 

APL = f ( eq A-2 

The friction factor, f, is estimated2 to be 96/Re. A 

value of 100/Re is used for this analysis. This is 

based on the flow along the tube bundles and the 

hydraulic diameter. Using the unobstructed 

velocity (VF) from the FLUENT predictions, 

Equation A-2 is reduced to Equation A-3. 

dPz= 770231 ., eq A-3 /VFe 

The viscous loss term used in the FLUENT porous model is set to 770 x 103 / M2 . T 1.173 

A.3.2 Inertial Loss r 

Inertial losses are accounted for at the lower 

orifice (nose), the entrance into the tube region (at 

lower tie plate), the seven grid spacers, and the 

upper tie plate final expansion. Values for these -*r- 0,4 

coefficients are determined from References 2 

and 33. At the low Reynolds (Re) numbers 

expected in the simulation, the loss coefficients lower i ?14 . 

vary strongly with Re. The uncertainty in these 

numbers is large. The net effect, however, is not N'-" 0 4F 

large, because the inertial losses are small e , " 

c o m p a re d to th e v is c o u s lo s s a t lo w R e . ..o.,.t 14..,3.t, l

Figure A-3. Bundle Geometry (meters)

A-3
Table A-2 lists the parameters used to represent the inertial losses within the bundles. The 'k' 

is the k factor obtained from the reference material. In some cases, interpolations or even 

extrapolations are made to obtain this value. A representative value of k is used when k varies 

strongly with Re. Vr represents the reference velocity description associated with the k factor. 

The terms kF and VF represent the friction factor and velocity based upon the FLUENT model 

flow area (100% open). The value C needed by the code is also given. Note that the effect of 

C in the code is not a series of abrupt pressure drops at physical rack features, but an 

equivalent total pressure drop spread over the entire rack length L. The net effect of the 

inertial loss coefficients is given in the last row of Table A-2. 

Table A-2. Inertial Loss Parameters 

k Vr based on k_ C = kF / L 

Lower orifice 11 flow in open channel box 19.11 4.45 

Entrance to tubes 1.2 flow along tubes 7.1 1.65 

Grid spacers (1 of 7) 3.2 flow along tubes 18.94 4.4 

Upper tie plate (w 4 flow along tubes 23.67 5.5 

expansion) j 

[Effect of lower orifice, entrance to tubes, 7 grid spacers, and upper tie plate J_42.5 / m 

A.3.3 Final Code Input Summary 

The FLUENT code input parameters which simulate the pressure drop through the rack and 

bundles are 770000 / m2 for viscous losses and 42.5 / m for inertial losses. These values apply 

in the vertical direction. Lateral resistance is infinite due to the solid rack walls. Using infinite 

loss coefficients is not recommended in the FLUENT code. The lateral resistance values used 

in the model are set two orders of magnitude higher than the vertical coefficients defined above. 

This increased resistance and the several vertical walls used in the CFD model align the flow in 

the vertical direction. 

A.4 COMPARISON WITH GE DATA 

A plot of pressure loss through an 8x8 GE bundle is obtained for a quick comparison with the 

results obtained above. The GE data represent a single-phase pressure drop through a BWR 

bundle as a function of flow rate. It is unclear if the GE results are measured or computed. 

The results are at low Re numbers similar to the Re numbers expected in this analysis. 

Although the geometry for the 8x8 bundle is different from the 9x9 geometry, these results are 

considered a good indication of the appropriateness of the loss coefficients determined above. 

The pressure drop for two different flow rates is computed from the coefficients given above 

and compared to the GE data in Figure A-4. The agreement is good. Units are omitted from 

the figure to protect potential proprietary information.

A-4
APPENDIX I-1: SINGLE PHASE BUNDLE 'aP IN STORAGE POOL 

10....... ____ D m o d el 

-mass flow rate 

Figure A-4. Comparison With Predicted Pressure Drop of GE Data 

(Units dropped to protect potential proprietary data) 

A.5 REFERENCES 

1. Fluent Inc., "Fluent 5 Users Guide," Vol 1, p 6-104, Canterra Resource Park, 10 

Cavendish Ct, Lebanon, NH 03766,1998. 

2. Todreas, N.E., Kazimi, M.S., Nuclear Systems /: Thermal Hydraulic Fundamentals, 

"Hemisphere Publishing Corporation, 1990. 

3. Idel'chik, 1. E., "Handbook of Hydraulic Resistance," AEC-TR-6630, 1960. 

A-5
APPENDIX B 

THERMAL PROPERTY DETERMINATION FOR POROUS 

REGION REPRESENTING BWR FUEL AND STORAGE RACKS 

B.1 BACKGROUND 

The computational fluid dynamics (CFD) analysis of spent fuel pool cooling documented in this 

report uses a porous medium to represent the racks and fuel, dispensing with the need to 

represent the geometric details of the rack and fuel structure. The porous model allows the 

user to input thermal properties for the racks and fuel to represent the thermal conduction and 

thermal inertia in this region. 

FLUENT's porous model uses an effective conductivity, as defined in equation B-11, to model 

the conduction in the porous media. 

keff = (Jkf +-(1- (D)ks eq B-1 

The fluid conductivity (kf) is augmented by the porous material conductivity (k.). The porosity of 

the medium is represented by 4). This appendix gives background information and 

assumptions used in determining ks and 0i. The fluid conductivity is obtained from a table for 

the thermal conductivity of air (used throughout the model). 

The thermal inertia term in the governing equations is also modified in the porous region to 

account for the properties of the solid material. Although thermal inertia only affects transient 

terms, its presence does stabilize the steady state solution. This appendix outlines the 

assumptions made to approximate the thermal inertia terms used in the FLUENT model. 

B.2 GEOMETRY 

The racks and fuel are considered a 

homogeneous nonisotropic structure whose ,2,8,in 

properties are defined by considering a single - q 134.0612 

rack cell. The cross-section of a single high- U 

density rack cell with a 9x9 BWR fuel bundle is Lx_ __m_ 9 

shown in Figure B-1. The dark regions on each (0.44 i oI 

side are boral inserts in the rack structure. The 0 0.566 in 

cell pitch is 0.154 m (6.06 in.). The fuel bundle braI • 0 , 

is contained in a channel box with an inner i 00 bo 1) = 13406 4 

box 00=1137.3653 dimension of 0.134 m (5.3 in.). There is a gap 00 

between the rack walls and the fuel channel box. 

For the purposes of this analysis, all fuel rods water rod 

are considered to be full length. The rod OD is 2141° t0, 0 0 

0.1118 m (.44 in.). The rod pitch is 0.1438 m 00 00 0. O 

(.566 in.). There are 74 fuel rods and 2 water 0000 00 00 , 

rods in the 9x9 cross-section. Table B-1 gives . 

other values used to describe the geometry for 

the purposes of determining the loss 149 .•ID I 

coefficients. - . 9 p.t1 

Figure B-1. Cross-Section of Bundle in Rack

B-1
The FLUENT code uses a single set of constant material properties and a porosity to account 

for the porous material. The rack and fuel structure in Figure B-1 is far more complex than the 

FLUENT model. For the purposes of this appendix, the rack and fuel are broken down into five 

regions according to material properties. The five regions are the rack, the air, the zirconium 

alloy, the fuel, and the boral inserts. Table B-1 lists the cross-sectional area (volume per meter 

length) and the thermal properties of the five regions. 

Table B-i. Material Properties and Area of the Five Material Regions 

Region Area (m2) p (kg/m3) Cp (J/kg-K) k (WIm2-K) 

rack 0.000939 7800 580 23 

air 0.013943 0.405 (600 °C) 1114 (600 0C) 0.061 (600 0C) 

zirconium 0.002274 6500 330 13 

fuel 0.006241 10970 247 3.6 

boral 0.000319 2500 1500 18 

Air, though included in the table, is not used in the analysis of thermal inertia. Its density is too 

low to have a measurable effect. Air is used in estimating the thermal conductivity but its effect 

is minimal. Properties for air are obtained from the table used in the FLUENT model for this 

case. Values at 600 0C are selected. This temperature is typical of the hottest fuel region in 

the racks. The rack is assumed to be stainless steel (Type 316). The steel properties come 

from Table 8-2 in Nuclear Systems 12. The fuel clad, water pipe, and channel box are assumed 

to consist of Zircaloy 2. Zircaloy properties are also obtained from Table 8-2 in Reference 2. 

The fuel (U0 2) properties are obtained from Table 8-1 in reference 2. The boral inserts are 

assumed to be boron carbide. The properties for this material are obtained from NUREG/ CR

61503. 

B.3 THERMAL INERTIA 

The thermal inertia term does not play a direct role in the final solution since it applies directly 

only to the transient terms. The final results are steady state. These terms do, however, help 

to stabilize the steady state results. Without this term, the fuel region temperature oscillations 

are very high due to the lack of thermal inertial in the model. For the purposes of estimating the 

thermal inertia, the transient term given by Equation B-2 is used. 

( 4Dpf Cpf +(1- 1)pS Cp,) dT / dt eq B-2 

Subscript f refers to the fluid (air) and subscript s refers to the solid. In this case, the solid is 

made up of several solids. The unknowns are 01, p., and Cps. These unknowns are 

determined as a composite of the materials given in Table B-I. Equation B-3 is used to find 

these values. 

(1- D) pCp, (., Cp, eq B-3

B-2
Equation B-3 states that the thermal inertial of the sought-after solid is equivalent to the sum of 

the thermal inertia of the solids in the rack and fuel structure. The summation is over the four 

materials defined in Table B-1 (air is excluded). The term (p refers to the volume fraction of the 

particular material. The term p Cp refers to the thermal inertia of each of the materials. The 

data are summarized in Table B-2. 

Table B-2. Thermal Inertia Terms 

Region % volume ((p) p Cp V (J 1K) 

(V is based on 1 m3 total Volume) 

rack 0.0396 179.1 x 103 

zirconium 0.0959 205.7 x 103 

fuel 0.2632 713.1 x 103 

boral 0.0135 50.4 x 103 

Summation 0.4121 1148 x 103 

The summation on the right side of equation B-3 is equal to 1148 x 10' J/K. At this point there 

is a single equation with three unknowns. The unknowns are the porosity, density, and specific 

heat on the left side of equation B-3. The properties of steel are selected (arbitrarily) and the 

porosity is computed to yield the correct thermal inertia. This results in a porosity of 0.746, 

which is rounded to 0.75, then applied to the model. 

Portions of the rack above and below the active fuel region do not have the same thermal 

inertia. As an approximation, the thermal inertia is computed assuming that the fuel and 

zirconium in Table B-2 are not present in these regions. Using the same material properties for 

the porous solid gives a porosity of 0.95 above and below the active fuel region. The thermal 

inertia parameters used in the model is given below. 

B.3.1 Input Parameters for Thermal Inertia 

"PS = 7800 kg/M3 

Cps = 580 J/kg-K 

S0 = 0.75 (within the fuel region) 

S0 = 0.95 (for the portion of the rack above and below the fuel) 

B.4 THERMAL CONDUCTION 

The thermal conduction plays a direct role in the final CFD model. This term directly influences 

a path of heat transfer away from the hottest fuel. A quick method is used to approximate this 

value, and sensitivity studies are done with the FLUENT model to quantify the effect of 

uncertainty in this input parameter. 

To estimate the thermal conduction, the 1 D parallel conduction resistance method is used. 

Equation B-4 defines the total resistance (RT) in terms of the component resistances (R,).

B-3
1 Y, 1 I eq B-4 

RT Ri 

Individual resistances are given by LI(k A). The term L is the distance in the conduction 

direction. The term k A is the thermal conductivity multiplied by the conduction area. Using the 

data from Table B-1 and assuming a length of 1 m gives a total resistance of 12.6. Using the 

total area of solids from Table B-1 (AT = 0.009733 M2 ), the effective conductivity can be 

computed from RT = 1/kAT. The effective conductivity is 8.12 for the solid region. This value 

also represents a simple volume-weighted thermal conductivity for the solid material regions. 

The value computed for the effective thermal conductivity must be modified to account for the 

actual porosity used in the model. The term {(1 -0) kJ} from equation B-1 is set equal to the 

conductivity found above and the porosity obtained by comparing the area fractions in Table B

1. The result is that the solid material thermal conductivity should be 13.4 for the porosity of 

0.75 set earlier. The porosity and/or thermal conductivity values are not significant by 

themselves. It is the combination {(1-0i) kJ} found in Equation B-1 that is significant. This 

appendix defines the porosity on the basis of the assumptions made in determining the thermal 

inertia terms. This porosity is carried through to the thermal conduction term since a region can 

only have one porosity. Even though the porosity value is not correct, the thermal conductivity 

is adjusted to yield the appropriate net effect. The thermal conduction parameters used in the 

model are given below. 

B.4.1 Summary of Input Parameters for Thermal Conductivity 

ks = 13.5 W/m-K (13.4 computed above; 13.5 used for model) 

0i' = 0.75 (in the fuel region) 

• 0' = 0.95 (for the portion of the rack above and below the fuel) 

B.5 REFERENCES 

1. Fluent Inc., "Fluent 5 Users Guide," Vol 1, p 6-107, Canterra Resource Park, 10 

Cavendish Ct, Lebanon, NH 03766, 1998. 

2. Todreas, N.E., Kazimi, M.S., Nuclear Systems /: Thermal Hydraulic Fundamentals, 

Hemisphere Publishing Corporation, 1990. 

3. SCDAP/RELAP5 Development Team, "SCDAP/RELAP5/MOD 3.2 Code Manual," 

NUREG/CR 6150, Vol. 4, Rev. 1, US Nuclear Regulatory Commission, Washington, DC 

20555, Jul1998.

B-4
NRC FORM 335 U.S. NUCLEAR REGULATORY COMMISSION 1. REPORT NUMBER 

(2-89) (Assigned by NRC, Add Vol., Supp., Rev., 

NRCM 1102, and Addendum Numbers, if any.) 

3201, 3202 BIBLIOGRAPHIC DATA SHEET 

(See instructions on the reverse) 

2. TITLE AND SUBTITLE NUREG-1726 

Predictions of Spent Fuel Heatup After a Complete Loss of Spent Fuel Pool Coolant 3. DATE REPORT PUBLISHED 

MONTH YEAR 

June 2000 

4. FIN OR GRANT NUMBER 

5. AUTHOR(S) 6. TYPE OF REPORT 

Christopher F. Boyd Technical 

7. PERIOD COVERED (Inclusive Dates) 

8. PERFORMING ORGANIZATION - NAME AND ADDRESS (If NRC, provide Division, Office or Region, U.S. Nuclear Regulatory Commission, end mailing address; if contractov, 

provide name end mailing address.) 

Division of Systems Analysis and Regulatory Effectiveness 

Office of Nuclear Regulatory Research 

U.S. Nuclear Regulatory Commission 

Washington, DC 20555-0001 

9. SPONSORING ORGANIZATION - NAME AND ADDRESS (If NRC, type 'Same as above"; if contractor, provide NRC Division, Office or Region, U.S. Nuclear Regulatory Commission, 

and mailing address.) 

Same as above

10. SUPPLEMENTARY NOTES 

11. ABSTRACT (200 words or less) 

This generic study was done to predict the peak fuel temperatures and flow patt erns observed after a complete loss of spent 

fuel pool coolant. The pool is assumed to be housed in a large containment bui Iding with an operable ventilation path. Fuel 

temperatures and natural circulation phenomena are predicted using computationa I fluid dynamics (CFD). These predictions 

identify the peak temperatures and the principal flow paths and features of the natural convective cooling phenomena for the 

fuel pool and containment building. The model is representative of a typical BWR pool. The pool is filled to capacity with 

4200 fuel bundles in high-density racking. Predictions are made with fuel load s representative of a fuel pool 2, 3, 4, and 6 

years after reactor shut down. Sensitivity studies are done on the ventilation rate, the outer wall heat transfer, the location of 

the fuel, fuel burnup, the flow resistance, and heat conduction within the rack s. The FLUENT CFD code predictions show that 

fuel temperatures remain below the temperature limits of 800 oC and 600 oC afte r 26 and 35 months, respectively. In 

considering the quantitative results of this study, the limitations and modelin g assumptions must be kept in mind.

12. KEY WORDS/DESCRIPTORS (List words or phrases that will assist researchers in locating the report.) 13. AVAILABILITY STATEMENT 

spent fuel pool unlimited 

complete loss of coolant 14. SECURITY CLASSIFICATION 

BWR (This Page) 

high density fuel racks unclassified 

fuel temperature (Tis Report) 

natural circulation flow unclassified 

decommissioning 15. NUMBER OF PAGES 

16. PRICE

NRC FORM 335 (2-89) This form was electronically produced by Elite Federal Forms, Inc.
NUREG-1726 PREDICTIONS OF SPENT FUEL HEATUP AFTER A COMPLETE LOSS OF SPENT FUEL POOL COOLANT JUNE 2000