Material Property Correlations: 

Comparisons between 

FRAPCON-3.4, FRAPTRAN 1.4, and

MATPRO

Office of Nuclear Regulatory Research 

NUREG/CR-7024 

PNNL-19417
Material Property Correlations:

Comparisons between

FRAPCON-3.4, FRAPTRAN 1.4, and

MATPRO

Manuscript Completed: August 2010 

Date Published: March 2011 

Prepared by

W.G. Luscher and K.J. Geelhood 

Pacific Northwest National Laboratory

P.O. Box 999

Richland, WA 99352

M.E. Flanagan, NRC Project Manager

NRC Job Code N6609

Office of Nuclear Regulatory Research 

NUREG/CR-7024 

PNNL-19417
iii 

Abstract 

The U.S. Nuclear Regulatory Commission uses the computer codes FRAPCON-3 and 

FRAPTRAN to model steady-state and transient fuel behavior, respectively, in regulatory 

analysis. In this document, material property correlations for oxide fuels and cladding materials 

are presented and discussed. Comparisons are made between the material property correlations 

used in the most recent versions of the codes, FRAPCON-3.4 and FRAPTRAN 1.4. 

Comparisons are also made with MATPRO, which is a compilation of material property 

correlations with an extensive history of use with various fuel performance and severe accident 

codes. In addition to model-to-model comparisons, model-to-data comparisons and source codes 

for the FRAPCON-3.4 and FRAPTRAN 1.4 correlations are provided. 

 
iv 
v 

Foreword 

The U.S. Nuclear Regulatory Commission uses the computer codes FRAPCON and FRAPTRAN 

to model steady-state and transient fuel behavior, respectively, in regulatory analysis. To 

effectively model fuel behavior, material property correlations must be used for a wide range of 

operating conditions (e.g., temperature and burnup). In this sense, a “material property” is a 

physical characteristic of the material whose quantitative value is necessary in the analysis 

process. 

The consolidated resource for “material properties” cited most often in the literature is MATPRO 

(as documented in Volume 4 of NUREG/CR-6150). MATPRO is a compilation of fuel and 

cladding material property correlations with an extensive history of use with fuel performance 

and severe accident codes. However, consistency among the material property correlations in 

FRAPCON-3, FRAPTRAN, and MATPRO has never been complete, as the material property 

correlations in the fuel performance codes have evolved to accommodate recent fuel design 

changes. 

The primary purpose of this report is to consolidate the current material property correlations 

used in FRAPCON and FRAPTRAN into a single document. In addition to describing the 

material property correlations used in the subroutines of latest versions of the codes, FRAPCON3.4 and FRAPTRAN 1.4, this report also provides a variety of comparisons to the material 

property correlations used in MATPRO. 

 
vi 
vii 

Contents 

Abstract .......................................................................................................................................... iii

Foreword .......................................................................................................................................... v

List of Figures ................................................................................................................................. ix

List of Tables ................................................................................................................................... xi

Executive Summary ..................................................................................................................... xiii

Abbreviations ................................................................................................................................. xv

1.0 Introduction ........................................................................................................................... 1.1

1.1 Relation to Other Reports .................................................................................................. 1.2

1.2 References .......................................................................................................................... 1.3

2.0 Oxide Fuel Properties ............................................................................................................ 2.1

2.1 Fuel Melting Temperature (PHYPRP) ............................................................................... 2.1

2.1.1 Model Development and Comparisons ...................................................................... 2.1

2.1.2 Applicability and Uncertainty .................................................................................... 2.3

2.1.3 References .................................................................................................................. 2.3

2.2 Fuel Specific Heat Capacity (FCP) and Fuel Enthalpy (FENTHL) ................................... 2.4

2.2.1 Model Development and Comparisons ...................................................................... 2.4

2.2.2 Applicability and Uncertainty .................................................................................... 2.7

2.2.3 References .................................................................................................................. 2.7

2.3 Fuel Thermal Conductivity (FTHCON)............................................................................. 2.8

2.3.1 Model Development and Comparisons ...................................................................... 2.8

2.3.2 Applicability and Uncertainty .................................................................................. 2.20

2.3.3 References ................................................................................................................ 2.21

2.4 Fuel Emissivity (FEMISS) ............................................................................................... 2.22

2.4.1 Model Development and Comparisons .................................................................... 2.22

2.4.2 Applicability and Uncertainty .................................................................................. 2.23

2.4.3 References ................................................................................................................ 2.23

2.5 Fuel Thermal Expansion (FTHEXP) ............................................................................... 2.24

2.5.1 Model Development and Comparisons .................................................................... 2.24

2.5.2 Applicability and Uncertainty .................................................................................. 2.26

2.5.3 References ................................................................................................................ 2.26

2.6 Fuel Densification (FUDENS) ......................................................................................... 2.27

2.6.1 Model Development and Comparisons .................................................................... 2.27

2.6.2 Applicability and Uncertainty .................................................................................. 2.30

2.6.3 References ................................................................................................................ 2.30

2.7 Fuel Swelling (FSWELL) ................................................................................................ 2.31

2.7.1 Model Development and Comparisons .................................................................... 2.31

2.7.2 Applicability and Uncertainty .................................................................................. 2.34

2.7.3 References ................................................................................................................ 2.34

3.0 Cladding Material Properties ................................................................................................ 3.1

3.1 Cladding Specific Heat (CCP) ........................................................................................... 3.1

3.1.1 Model Development and Comparisons ...................................................................... 3.1

3.1.2 Applicability and Uncertainty .................................................................................... 3.3

3.1.3 References .................................................................................................................. 3.3

3.2 Cladding Thermal Conductivity (CTHCON) ..................................................................... 3.3

3.2.1 Model Development and Comparisons ...................................................................... 3.4

3.2.2 Applicability and Uncertainty .................................................................................... 3.6

3.2.3 References .................................................................................................................. 3.6

3.3 Cladding Oxide Thermal Conductivity (ZOTCON) .......................................................... 3.7
viii 

3.3.1 Model Development and Comparisons ...................................................................... 3.8

3.3.2 Applicability and Uncertainty .................................................................................. 3.11

3.3.3 References ................................................................................................................ 3.11

3.4 Cladding Surface Emissivity (ZOEMIS) ......................................................................... 3.11

3.4.1 Model Development and Comparisons .................................................................... 3.12

3.4.2 Applicability and Uncertainty .................................................................................. 3.14

3.4.3 References ................................................................................................................ 3.14

3.5 Cladding Thermal Expansion (CTHEXP) ....................................................................... 3.14

3.5.1 Model Development and Comparisons .................................................................... 3.14

3.5.2 Applicability and Uncertainty .................................................................................. 3.17

3.5.3 References ................................................................................................................ 3.18

3.6 Cladding Elastic Modulus (CELMOD) and Shear Modulus (CSHEAR) ........................ 3.18

3.6.1 Model Development and Comparisons .................................................................... 3.19

3.6.2 Applicability and Uncertainty .................................................................................. 3.24

3.6.3 References ................................................................................................................ 3.24

3.7 Cladding Axial Growth (CAGROW) .............................................................................. 3.25

3.7.1 Model Development and Comparisons .................................................................... 3.25

3.7.2 Applicability and Uncertainty .................................................................................. 3.29

3.7.3 References ................................................................................................................ 3.29

3.8 Creep Rate (CREPR) ....................................................................................................... 3.30

3.8.1 Model Development and Comparisons .................................................................... 3.30

3.8.2 Applicability and Uncertainty .................................................................................. 3.33

3.8.3 References ................................................................................................................ 3.34

3.9 Cladding Meyer Hardness (CMHARD)........................................................................... 3.34

3.9.1 Model Development and Comparisons .................................................................... 3.35

3.9.2 Applicability and Uncertainty .................................................................................. 3.36

3.9.3 References ................................................................................................................ 3.36

4.0 Gas Material Properties ......................................................................................................... 4.1

4.1 Gas Conductivity (GTHCON) ........................................................................................... 4.1

4.1.1 Model Development and Comparisons ...................................................................... 4.1

4.1.2 Applicability and Uncertainty .................................................................................... 4.7

4.1.3 References .................................................................................................................. 4.8

Appendix A Subroutine Source Codes ......................................................................................... A.1
ix 

List of Figures 

Figure 2.1-1 A comparison between solidus temperatures calculated by MATPRO using the 

FHYPRP subroutine and FRAPCON-3/FRAPTRAN (FRAP) using the 

PHYPRP subroutine. Results for burnup varying from 0 to 62 GWd/MTU and 

Pu contents of 0, 3, and 7 wt% are presented. .................................................. 2.2

Figure 2.1-2 A comparison between liquidus temperatures calculated by MATPRO using 

the FHYPRP subroutine and FRAPCON-3/FRAPTRAN (FRAP) using the 

PHYPRP subroutine. Results for burnup varying from 0 to 62 GWd/MTU and 

Pu contents of 0, 3, and 7 wt% are presented. .................................................. 2.3

Figure 2.2-1 Specific heat capacity of UO2 from three experimenters compared with the 

FRAPCON-3/FRAPTRAN/MATPRO (FRAP) correlation (solid line) for UO2.2.6

Figure 2.2-2 Specific heat capacity of (U0.8Pu0.2)O2 from three experimenters compared with 

the FRAPCON-3/FRAPTRAN/MATPRO (FRAP) correlations (solid line) for 

MOX. ............................................................................................................... 2.7

Figure 2.3.3-1 MATPRO and modified NFI thermal conductivity models for UO2 as a 

function of temperature. Unlike the MATPRO model, the modified NFI model 

can account for burnup and is presented at burnup levels of 0 and 

30 GWd/MTU for comparison. ...................................................................... 2.14

Figure 2.3.3-2 MATPRO model-to-data comparison for thermal conductivity of unirradiated 

UO2. ................................................................................................................ 2.15

Figure 2.3.3-3 Modified NFI model-to-data comparison for thermal conductivity of 

unirradiated UO2. ........................................................................................... 2.15

Figure 2.3.3-4 MATPRO model-to-data comparison for thermal conductivity of irradiated 

UO2. ................................................................................................................ 2.16

Figure 2.3.3-5 Modified NFI model-to-data comparison for thermal conductivity of irradiated 

UO2. ................................................................................................................ 2.17

Figure 2.3.3-6 MATPRO and Duriez/modified NFI thermal conductivity models as a function 

of temperature for MOX fuel with 7 wt% plutonium content. Unlike the 

MATPRO model, the Duriez/modified NFI model can account for burnup and 

is presented at burnup levels of 0 and 30 GWd/MTU for comparison. ......... 2.18

Figure 2.3.3-7 MATPRO model-to-data comparison for thermal conductivity of MOX 

modeled with 7 wt% of plutonium. ................................................................ 2.19

Figure 2.3.3-8 Duriez/modified NFI model-to-data comparison for thermal conductivity of 

MOX with modeled with 7 wt% of plutonium. ............................................. 2.19

Figure 2.4-1 Data to model comparison for FEMISS. ........................................................ 2.23

Figure 2.5-1 Model-to-data comparison for UO2 correlation used in MATPRO and 

FRAPCON-3/FRAPTRAN with updated constants. ..................................... 2.25

Figure 2.5-2 Model-to-data comparison for PuO2 correlation used in MATPRO and 

FRAPCON-3/FRAPTRAN. ........................................................................... 2.26

Figure 2.6-1 Comparison of densification predictions generated by the FUDENS model and 

experimental densification measurements determined in reactor. ................. 2.30

Figure 2.7-1 Data comparison between MATPRO and FRAPCON-3.4 versions of 

FSWELL. ....................................................................................................... 2.33

Figure 2.7-2 Comparison between swelling rates determined for various instrumented fuel 

assemblies tests at Halden and swelling rates predicted by FRAPCON-3.4. . 2.34

Figure 3.2-1 CTHCON model and data used for development. ........................................... 3.5

Figure 3.2-2 CTHCON model and data collected since development. ................................. 3.6
x 

Figure 3.3-1 Oxide thermal conductivity measurements as a function of temperature 

compared with the MATPRO and FRAPCON-3.4/FRAPTRAN 1.4 (FRAP) 

cladding oxide thermal conductivity models. ................................................ 3.10

Figure 3.4-1 Model-to-data comparison for cladding oxide emissivity. ............................. 3.13

Figure 3.5-1 Model-to-data comparison for FRAPCON-3.4/FRAPTRAN 1.4 (FRAP) and 

MATPRO axial thermal expansion. Data from cladding tubes..................... 3.17

Figure 3.5-2 Model-to-data comparison for FRAPCON-3.4/FRAPTRAN 1.4 (FRAP) and 

MATPRO diametral thermal expansion. Data from cladding tubes. ............ 3.17

Figure 3.6-1 Predicted vs. measured values for Young’s modulus for various data sets. .. 3.23

Figure 3.7-5 Model-to-data comparison for SRA Zircaloy-4 (PWR). The standard error of 

L/L (%) is 0.11. ............................................................................................ 3.27

Figure 3.7-6 Model-to-data comparison for RXA Zircaloy-2 (BWR). The standard error of 

L/L (%) is 0.07. ............................................................................................ 3.27

Figure 3.7-7 Model-to-data comparison for M5. M5 data was collected by Gilbon et al. 

(2000). The standard error of L/L (%) is 0.05. ........................................... 3.28

Figure 3.7-8 Model-to-data comparison for ZIRLO. Vandellos data was collected by Irisa 

(2000) while North Anna and BR-3 data was collected by Sabol (1994). The 

standard error of L/L (%) is 0.05. ................................................................ 3.28

Figure 3.8-1 Predicted vs. measured hoop strain for SRA Zircaloy model. ....................... 3.33

Figure 3.8-2 Predicted vs. measured hoop strain for RXA Zircaloy model. ...................... 3.33

Figure 3.9-1 CMHARD correlation as a function of temperature with data from Peggs and 

Godin. ............................................................................................................. 3.36

Figure 4.1-1 Model-to-data comparison for He thermal conductivity. ................................. 4.3

Figure 4.1-2 Model-to-data comparison for Ar thermal conductivity. ................................. 4.4

Figure 4.1-3 Model-to-data comparison for Kr thermal conductivity. ................................. 4.5

Figure 4.1-4 Model-to-data comparison for Xe thermal conductivity. ................................. 4.6

Figure 4.1-5 Model-to-data comparison for H2 thermal conductivity. ................................. 4.6

Figure 4.1-6 Model-to-data comparison for N2 thermal conductivity. ................................. 4.7

Figure 4.1-7 Model-to-data comparison for steam thermal conductivity. ............................ 4.7
xi 

List of Tables 

Table 1.1 Roadmap to documentation of models and properties in NRC fuel performance 

codes, FRAPCON-3.4 and FRAPTRAN 1.4. .................................................. 1.2

Table 2.2-1 Constants used in UO2, Gd2O3, and PuO2 heat capacity and enthalpy 

correlations. ...................................................................................................... 2.5

Table 2.5-1 Parameters used in UO2 and PuO2 solid-phase thermal expansion correlations.2.24

Table 3.1-1 Specific heat capacity database for Zircaloy-2, Zircaloy-4, ZIRLO, and M5. . 3.2

Table 3.1-2 Specific heat capacity database for ZrNb-1 alloys. .......................................... 3.2

Table 3.8-1 Parameters for FRAPCON-3.4 creep equation for SRA and RXA cladding. 3.31

Table 4.1-1 Constants used in gas thermal conductivity correlations in MATPRO. ........... 4.1

Table 4.1-2 Constants used in gas thermal conductivity correlations in FRAP. ................. 4.2

Table 4.1-3 Standard deviation of gas thermal conductivity correlations. .......................... 4.8
xii 
xiii 

Executive Summary 

The U.S. Nuclear Regulatory Commission uses the computer codes FRAPCON-3 and 

FRAPTRAN to model steady-state and transient fuel behavior, respectively, in regulatory 

analysis. In this document, material property correlations for oxide fuels and cladding materials 

are presented and discussed. Comparisons are made between the material property correlations 

used in the most recent versions of the codes, FRAPCON-3.4 and FRAPTRAN 1.4. 

Comparisons are also made with MATPRO, which is a compilation of material property 

correlations with an extensive history of use with various fuel performance and severe accident 

codes. In addition to model-to-model comparisons, model-to-data comparisons and source codes 

for the FRAPCON-3.4 and FRAPTRAN 1.4 correlations are also provided. 

This is one of three documents that describe the FRAPCON-3.4 and FRAPTRAN 1.4 fuel 

performance codes. The FRAPCON-3.4 code structure and behavioral models are described in 

the FRAPCON-3.4 code description document (Geelhood et al., 2010a). The FRAPTRAN 1.4 

code structure and behavioral models are described in the FRAPTRAN 1.4 code description 

document (Geelhood et al., 2010b). Additional descriptions of mechanical models and model-todata comparisons for mechanical property correlations are provided in the mechanical properties 

report (Geelhood, 2008). 
xiv 
xv 

Abbreviations 

°C degrees Celsius 

°F degrees Fahrenheit 

Ar argon 

BWR boiling-water reactor 

cm2

 square centimeter(s) 

crud Chalk River Unidentified Deposit (generic term for various residues deposited on 

fuel rod surfaces, originally coined by Atomic Energy of Canada, Ltd. to describe 

deposits observed on fuel from the test reactor at Chalk River) 

Gd2O3 gadolinia 

GPa gigapascal(s) 

GWd/MTU gigawatt-days per metric ton of uranium 

H hydrogen 

He helium 

J joule(s) 

K kelvin 

kg kilogram(s) 

Kr krypton 

LWR light-water reactor 

m meter(s) 

m2

 square meter(s) 

m3

 cubic meter(s) 

MeV megaelectron volt(s) 

mm millimeter(s) 

MOX mixed oxide, (U, Pu)O2

MPa megapascal(s) 

MWd/kgU megawatt-days per kilogram of uranium 

MWd/MTU megawatt-days per metric ton of uranium 

MWs/kgU megawatt-seconds per kilogram of uranium 

n neutron(s) 

N nitrogen 

Nb niobium

NFI Nuclear Fuels Industries 

NRC U.S. Nuclear Regulatory Commission 

O/M oxygen-to-metal 

Pa pascal(s) 

PNNL Pacific Northwest National Laboratory 

Pu plutonium 

PuO2 plutonium dioxide 

PWR pressurized-water reactor 

RXA recrystallized annealed 

s second(s) 

SRA stress relief annealed 

TD theoretical density 

tHM tonne(s) of heavy metal 

UO2 uranium dioxide 

UO2-Gd2O3 urania-gadolinia 

W watt(s) 

wt% percent by weight 
xvi 

Xe xenon 

ZrNb-1 Russian zirconium alloy with 1 wt% of niobium 

ZrO2 zirconium dioxide 

m micrometer(s) 
1.1 

1.0 Introduction 

The U.S. Nuclear Regulatory Commission (NRC) uses the computer codes FRAPCON-3 and 

FRAPTRAN to model steady-state and transient fuel behavior, respectively, in regulatory 

analysis. To effectively model fuel behavior, material property correlations must be used for a 

wide range of operating conditions (e.g., temperature and burnup). In this sense, a “material 

property” is a physical characteristic of the material whose quantitative value is necessary in the 

analysis process. Further, the property may be used to compare the benefits of one material with 

those of another. Generally speaking, the material properties of interest in regulatory analysis of 

nuclear fuel behavior are mechanical or thermodynamic. 

The issue of what is and is not a “material property” will never be universally resolved. In this 

report, properties such as thermal conductivity are included. Other characteristics of the material 

(e.g., fission gas release) are considered “models” rather than properties and are discussed 

elsewhere. The material property correlations for the FRAPCON-3 and FRAPTRAN computer 

codes were partially documented in NUREG/CR-6534 and NUREG/CR-6739, respectively 

(Lanning et al., 2005; Cunningham et al., 2001) as well as in NUREG/CR-6150 (Siefken et al., 

2001). Some of these have been modified or updated since the original code documentation was 

published. 

The primary purpose of this report is to consolidate the current material property correlations 

used in FRAPCON-3 and FRAPTRAN into one document. Material property correlations for 

oxide fuels, including uranium dioxide (UO2) and mixed oxide (MOX) fuels, are described in 

Section 2. Throughout this document, the term MOX is used to describe fuels that are blends of 

uranium and plutonium oxides, (U,Pu)O2. The properties for UO2 with other additives 

(e.g., gadolinia) are also discussed. Material property correlations for cladding materials and 

gases are described in Sections 3 and 4, respectively. 

In addition to describing the material property correlations used in the subroutines of 

FRAPCON-3 and FRAPTRAN, this report provides a variety of comparisons between material 

property correlations and data. Although they are frequently identical, comparisons are made 

between the material property correlations used in the FRAPCON-3 and FRAPTRAN codes. 

Comparisons are also made between the material property correlations used in MATPRO, a 

compilation of fuel and cladding material property correlations with an extensive history of use 

with various fuel performance and severe accident codes. For several reasons, consistency 

among the material property correlations in FRAPCON-3, FRAPTRAN, and MATPRO has never 

been complete. However, the current versions of FRAPCON-3 and FRAPTRAN use a relatively 

consistent set of correlations for the properties that are used by both codes. The material property 

correlations in the most recent version of MATPRO are documented in Volume 4 of 

NUREG/CR-6150 (Siefken et al., 2001). In addition to comparison of the various correlations, 

correlation-to-data comparisons are made with FRAPCON-3, FRAPTRAN, and MATPRO 

properties. 

All comparisons in this report are based on the material property correlations used in the most 

recent version of the FRAPCON-3 and FRAPTRAN codes, FRAPCON-3.4 and FRAPTRAN 1.4. 

The source code for each material property correlation discussed is provided for FRAPCON-3.4 

and FRAPTRAN 1.4 (see Appendix A) as well as a range of applicability and an estimate of 

uncertainty, where possible. 
1.2 

1.1 Relation to Other Reports 

The full documentation of the steady-state and transient fuel performance codes is described in 

three documents. The basic fuel, cladding, and gas material properties used in FRAPCON-3.4 

and FRAPTRAN 1.4 are described in the material properties handbook (this report). The 

FRAPCON-3.4 code structure and behavioral models are described in the FRAPCON-3.4 code 

description document (Geelhood et al., 2010a). The FRAPTRAN 1.4 code structure and 

behavioral models are described in the FRAPTRAN 1.4 code description document (Geelhood et 

al., 2010b). Additional descriptions of mechanical models and model-to-data comparisons for 

mechanical property correlations are provided in the mechanical properties report (Geelhood, 

2008). 

Table 1.1 shows where each specific material property and model used in the NRC fuel 

performance codes are documented. 

Table 1.1 Roadmap to documentation of models and properties in NRC fuel 

performance codes, FRAPCON-3.4 and FRAPTRAN 1.4. 

Model/Property FRAPCON-3.4 FRAPTRAN 1.4 

Fuel thermal conductivity Material properties handbook Material properties handbook 

Fuel thermal expansion Material properties handbook Material properties handbook 

Fuel melting temperature Material properties handbook Material properties handbook 

Fuel specific heat Material properties handbook Material properties handbook 

Fuel enthalpy Material properties handbook Material properties handbook 

Fuel emissivity Material properties handbook Material properties handbook 

Fuel densification Material properties handbook NA 

Fuel swelling Material properties handbook NA 

Fission gas release FRAPCON-3 code description FRAPTRAN code description 

Fuel relocation FRAPCON-3 code description FRAPTRAN code description 

Fuel grain growth FRAPCON-3 code description NA 

High burnup rim model FRAPCON-3 code description NA 

Nitrogen release FRAPCON-3 code description NA 

Helium release FRAPCON-3 code description NA 

Radial power profile FRAPCON-3 code description NA (input parameter) 

Stored energy FRAPCON-3 code description FRAPTRAN code description 

Decay heat model NA FRAPTRAN code description 

Fuel and cladding temperature 

solution 

FRAPCON-3 code description FRAPTRAN code description 

Cladding thermal conductivity Material properties handbook Material properties handbook 

Cladding thermal expansion Material properties handbook Material properties handbook 

Cladding elastic modulus Material properties handbook Material properties handbook 

Cladding creep model Material properties handbook NA 

Cladding specific heat Material properties handbook Material properties handbook 
1.3 

Model/Property FRAPCON-3.4 FRAPTRAN 1.4 

Cladding emissivity Material properties handbook Material properties handbook 

Cladding axial growth Material properties handbook NA 

Cladding Meyer hardness Material properties handbook Material properties handbook 

Cladding annealing FRAPCON-3 code description FRAPTRAN code description 

Cladding yield stress, ultimate 

stress, and plastic deformation 

FRAPCON-3 code description FRAPTRAN code description 

Cladding failure criteria NA FRAPTRAN code description 

Cladding waterside corrosion FRAPCON-3 code description NA (input parameter) 

Cladding hydrogen pickup FRAPCON-3 code description NA (input parameter) 

Cladding high temperature 

oxidation 

NA FRAPTRAN code description 

Cladding ballooning model NA FRAPTRAN code description 

Cladding mechanical deformation FRAPCON-3 code description FRAPTRAN code description 

Oxide thermal conductivity Material properties handbook Material properties handbook 

Crud thermal conductivity FRAPCON-3 code description NA 

Gas conductivity Material properties handbook Material properties handbook 

Gap conductance FRAPCON-3 code description FRAPTRAN code description 

Plenum gas temperature FRAPCON-3 code description FRAPTRAN code description 

Rod internal pressure FRAPCON-3 code description FRAPTRAN code description 

Coolant temperature and heat 

transfer coefficients 

FRAPCON-3 code description FRAPTRAN code description 

Optional models and properties not developed at PNNL 

VVER fuel and cladding models NA NUREG/IA-0164 

(Shestopalov et al., 1999) 

Cladding FEA model VTT-R-11337-06 

(Knuttilla, 2006) 

VTT-R-11337-06 

(Knuttilla, 2006) 

NA = not applicable 

VVER = water-cooled, water-moderated energy reactor 

FEA = finite element analysis 

1.2 References 

Lanning, D.D., C.E. Beyer, and K.J. Geelhood. 2005. FRAPCON-3 Updates, Including MixedOxide Fuel Properties, NUREG/CR-6534, Vol. 4, PNNL-11513, Pacific Northwest National 

Laboratory, Richland, WA. 

Cunningham, M.E., C.E. Beyer, F.E. Panisko, P.G. Medvedev, G.A. Berna, and H.H. Scott. 

2001. FRAPTRAN: Integral Assessment, NUREG/CR-6739, Vol. 2, PNNL-13576, Pacific 

Northwest National Laboratory, Richland, WA. 
1.4 

Siefken, L.J., E.W. Coryell, E.A. Harvego, and J.K. Hohorst. 2001. SCDAP/RELAP5/MOD 3.3 

Code Manual MATPRO-A Library of Materials Properties for Light-Water-Reactor-Accident 

Analysis, NUREG/CR-6150, Vol. 4, Rev. 2, INEL-96/0422, Idaho National Engineering and 

Environmental Laboratory, Idaho Falls, ID. 

Geelhood, K.J., W.G. Luscher, and C.E. Beyer. 2010a. FRAPCON-3.4: A Computer Code for 

the Calculation of Steady-State, Thermal-Mechanical Behavior of Oxide Fuel Rods for High 

Burnup, NUREG/CR-7022, Vol. 1, PNNL-19418, Vol. 1, Pacific Northwest National Laboratory, 

Richland, WA. 

Geelhood, K.J., W.G. Luscher, C.E. Beyer, and J.M. Cuta. 2010b. FRAPTRAN 1.4: A Computer 

Code for the Transient Analysis of Oxide Fuel Rods, NUREG/CR-7023, Vol. 1, PNNL-19400, 

Vol. 1, Pacific Northwest National Laboratory, Richland, WA. 

Geelhood, K.J., W.G. Luscher, and C.E. Beyer. 2008. PNNL Stress/Strain Correlation for 

Zircaloy, PNNL-17700, Pacific Northwest National Laboratory, Richland, WA. 

Shestopalov, A., K. Lioutov, L. Yegorova, G. Abyshov, and K. Mikitiouk. 1999. Modification of 

USNRC’s FRAP-T6 Fuel Rod Transient Code for High Burnup VVER Fuel, NUREG/IA-0164, 

U.S. Nuclear Regulatory Commission, Washington, D.C. 

Knuutila, A. 2006. Improvements on FRAPCON-3/FRAPTRAN Mechanical Modeling, VTT-R11337-06, VTT, Finland.
2.1 

2.0 Oxide Fuel Properties 

Material property correlations for UO2 and MOX fuel pellets are described in the following 

subsections. These correlations include the melting temperature, specific heat capacity and 

enthalpy, thermal conductivity, emissivity, thermal expansion, densification, and swelling of the 

fuel pellets. When indicated, some correlations account for gadolinia (Gd2O3) additions in UO2

fuel pellets. The subroutine containing the material correlation in the FRAPCON-3/FRAPTRAN 

codes is given in parentheses in the heading of each subsection. 

2.1 Fuel Melting Temperature (PHYPRP) 

The subroutine PHYPRP is used to calculate the temperature of the first appearance of liquid 

phase (solidus) and the temperature of the melting point of the last solid phase (liquidus) of UO2

and MOX. These temperatures are calculated as a function of burnup and plutonia content. The 

same subroutine is used in both FRAPCON-3.4 and FRAPTRAN 1.4 and exhibits smaller burnup 

dependence than the corresponding MATPRO correlation, FHYPRP. 

2.1.1 Model Development and Comparisons 

The equations used to calculate the UO2 and MOX melting points were derived by using 

3113.15K as the melting temperature of urania, which was determined experimentally by 

Brassfield (1968), and a least squares fit to parabolic equations for the solidus and liquidus 

boundaries from the Lyon and Baily (1967) phase diagram for the stoichiometric (U, Pu)O2

mixed oxide. These correlations are presented in Equations 2.1-1 through 2.1-5. 

(2.1-1) 

(2.1-2) 

For c = 0 

(2.1-3) 

For c > 0 

(2.1-4) 

(2.1-5) 

where 

c = Pu content (wt%) 

sldus (c) = solidus temperature as a function of Pu content (°C) 

liqdus (c) = liquidus temperature as a function of Pu content (°C) 

ftmelt = fuel melting temperature (K) 

FBu = burnup (MWd/MTU) 

fdelta = temperature range between solidus and liquidus (K) 

23 c10x468390.7c41395.52840)c(sldus   

22 c10x448518.1c21860.32840)c(liqdus   

ftmelt   10000/FBu515.3113

fdelta    10000/FBu5)c(sldus)c(liqdus

ftmelt    10000/FBu515.273)c(sldus
2.2 

The correlations defined above for the subroutine PHYPRP are almost identical to those used in 

the subroutine FHYPRP documented in MATPRO. However, the burnup dependence of the 

PHYPRP subroutine differs from the FHYPRP subroutine. A burnup dependence of 

3.2 K/GWd/MTU is used in FHYPRP, while a burnup dependence of 0.5 K/GWd/MTU is used in 

PHYPRP. This change was suggested by Popov et al. (2000) and is based on data from Adamson 

et al. (1985) and Komatsu et al. (1988). 

The burnup-dependent term in Equations 2.1-3 through 2.1-5 reflects this change in burnup 

dependence. Also, the term 3.21860 c in Equation 2.1-2 is slightly different from the term used 

in FHYPRP. The coefficient was changed from 3.21660 in FHYPRP to 3.21860 in PHYPRP. 

However, this change does not appreciably affect the liquidus temperature. 

Comparisons between solidus and liquidus temperature calculated by PHYPRP and FHYPRP are 

presented in Figures 2.1-1 and 2.1-2. Both subroutines predict solidus and liquidus temperatures 

that decrease with increasing burnup and plutonium concentration. However, these figures reveal 

a smaller decline in solidus and liquidus temperatures for the FRAPCON-3 and FRAPTRAN 

subroutine (PHYPRP) than the MATPRO subroutine (FHYPRP). This is a result of reducing the 

burnup dependence from 3.2 to 0.5 K/GWd/MTU between the FHYPRP and PHYPRP 

subroutines, respectively. 

Solidus Temperatures

2825

2875

2925

2975

3025

3075

3125

0 10 20 30 40 50 60

Burnup (GWd/MTU)

Temperature (K)

MATPRO 0 wt% MATPRO 3 wt% MATPRO 7 wt%

FRAP 0 wt% FRAP 3 wt% FRAP 7 wt%

Figure 2.1-1 A comparison between solidus temperatures calculated by MATPRO using 

the FHYPRP subroutine and FRAPCON-3/FRAPTRAN (FRAP) using the 

PHYPRP subroutine. Results for burnup varying from 0 to 62 GWd/MTU 

and Pu contents of 0, 3, and 7 wt% are presented. 
2.3 

Liquidus Temperatures

2825

2875

2925

2975

3025

3075

3125

0 10 20 30 40 50 60

Burnup (GWd/MTU)

Temperature (K)

MATPRO 0 wt% MATPRO 3 wt% MATPRO 7 wt%

FRAP 0 wt% FRAP 3 wt% FRAP 7 wt%

Figure 2.1-2 A comparison between liquidus temperatures calculated by MATPRO using 

the FHYPRP subroutine and FRAPCON-3/FRAPTRAN (FRAP) using the 

PHYPRP subroutine. Results for burnup varying from 0 to 62 GWd/MTU 

and Pu contents of 0, 3, and 7 wt% are presented. 

2.1.2 Applicability and Uncertainty 

The correlation used in PHYPRP is valid through the melting temperature of UO2 and MOX for 

burnups up to 62 GWd/MTU. Currently, there are no models for urania-gadolinia (UO2-Gd2O3) 

fuel melting. However, there are plans to develop a correlation in a future revision of the code. 

Due to the experimental data and data interpretations on which these temperature correlations are 

based, it is difficult to express a quantitative measure of uncertainty. A conservative estimate of 

uncertainty based on engineering judgment is ± 50K. 

2.1.3 References 

Brassfield, H.C., et al. 1968. Recommended Property and Reactor Kinetics Data for Use in 

Evaluating a Light-Water-Coolant Reactor Loss-of-Coolant Incident Involving Zircaloy-4 or 304-

SS-Clad UO2, GEMP-482, General Electric Company, Cincinnati, OH. 

Lyon, W.F., and W.E. Baily. 1967. “The Solid Liquid Phase Diagram for the UO2-PuO2

System,” Journal of Nuclear Materials, 22:332. 

Popov, S.G., J.J. Carbajo, V.K. Ivanov, and G.L. Yoder. 2000. Thermophysical Properties of 

MOX and UO2 Fuel Including the Effects for Irradiation, ORNL/TM 2000/351, Oak Ridge 

National Laboratory, Oak Ridge, TN. 

Adamson, M.G., E.A. Aitken, and R.W. Caputi. 1985. “Experimental and thermodynamic 

evaluation of the melting behavior of irradiated oxide fuels.” Journal of Nuclear Materials, 

130:349-365. 
2.4 

Komatsu, J., T. Tachibana, and K. Konashi. 1988. “The Melting Temperature of Irradiated 

Oxide Fuel,” Journal of Nuclear Materials, 154:38-44. 

2.2 Fuel Specific Heat Capacity (FCP) and Fuel Enthalpy (FENTHL) 

The subroutines FCP and FENTHL are used to calculate the specific heat capacity and enthalpy 

of the fuel pellet, respectively. The specific heat capacity and enthalpy of nuclear fuel are 

modeled empirically as functions of four parameters: temperature, composition, molten fraction, 

and oxygen-to-metal (O/M) ratio. The same subroutine is used in FRAPCON-3.4, FRAPTRAN 

1.4, and MATPRO. 

2.2.1 Model Development and Comparisons 

Equations for the specific heat capacity and enthalpy of solid UO2 and plutonium dioxide (PuO2) 

are assumed to have the same form, but with different constants. The basic relationships are 

presented in Equations 2.2-1 and 2.2-2, respectively. 

(2.2-1) 

(2.2-2) 

where 

FCP = specific heat capacity (J/kg*K) 

FENTHL = fuel enthalpy (J/kg) 

T = temperature (K) 

Y = oxygen-to-metal ratio 

R = universal gas constant = 8.3143 (J/mol*K) 

 = the Einstein temperature (K) 

ED = activation energy for Frenkel defects (J/mol) 

and the constants are given in Table 2.2-1. 

  )RT/E(exp RT2

EYK

TK

1)T/(expT

)T/(expK FCP 2 D

D3

2 2 2

2

1   

 

 )RT/E(expK 2

Y

2

TK

1)T/(exp

K FENTHL 3 D

2

1 2   

 
2.5 

Table 2.2-1 Constants used in UO2, Gd2O3, and PuO2 heat capacity and enthalpy 

correlations. 

Constant UO2 PuO2 GdO2

* Units 

K1 296.7 347.4 315.86 J/kg*K 

K2 2.43 x 10-2 3.95 x 10-4 4.044 x 10-2 J/kg*K2

K3 8.745 x 107 3.860 x 107 0.0 J/kg 

 535.285 571.000 348.0 K 

ED 1.577 x 105 1.967 x 105 0.0 J/mol 

* Gadolinia additions are only accounted for in FRAPCON-3.4. 

The constants in Table 2.2-1 were determined by Kerrisk and Clifton (1972) for UO2 and Kruger 

and Savage (1968) for PuO2. The specific heat capacity of UO2 in the liquid state (see Equation 

2.2-3) was determined by Leibotwitz (1971) and assumed to be valid for PuO2 in the liquid state. 

 FCPLiquid = 503 J/kg*K (2.2-3) 

Inspection of Equations 2.2-1 and 2.2-2 reveals that fuel enthalpy, or stored energy, is calculated 

by integrating the fuel specific heat capacity correlation from 0 to T in kelvin. Stored energy is 

important in reactor transient analysis because the severity of the transient is significantly 

affected by the initial stored energy of the fuel. Since UO2 and PuO2 are the principal oxides in 

light-water reactor (LWR) fuels, they are the constituents considered. 

The correlations for specific heat and enthalpy in both UO2 and MOX fuel pellets are identical in 

FRAPCON-3.4, FRAPTRAN 1.4, and MATPRO. However, FRAPCON-3.4 can accommodate 

gadolinia additions whereas FRAPTRAN cannot. For a mixture of UO2, Gd2O3, and PuO2, the 

specific heat capacity of the solid is determined by combining the contribution from each 

constituent in proportion to its weight fraction. When the material is partially molten, the heat 

capacity is determined similarly with a weighted sum of the solid and molten fractions. 

Since the specific heat correlation is only valid above a fuel temperature of about 300K, the fuel 

enthalpy correlation is not valid below a temperature of about 300K. In addition, it is necessary 

to calculate fuel enthalpy with respect to a reference temperature  300K. Thus, the fuel enthalpy 

at any desired temperature, T, is calculated by evaluating Equation 2.2-2 at T and a reference 

temperature, Tref, of 300K and taking the difference [FENTHL(T) – FENTHL(Tref)]. For 

temperatures greater than 2K below melting, the molten fraction and heat of fusion are used to 

interpolate between the enthalpy of unmelted fuel and just-melted fuel at the melting temperature. 

Specific heat correlations are compared with independent experimental data in Figure 2.2-1 and 

Figure 2.2-2 for UO2 and PuO2, respectively. Comparison between the UO2 model and 

experimental data from Hein (1968), Leibowitz (1969), and Gronvold (1970) in Figure 2.2-1 

reveals agreement between predicted and measured values up to high temperatures (a few 

hundred kelvin below melting). At this point the data begin to fall lower than the model. This is 

attributed to partial melting due to a non-uniform temperature distribution within the sample. 

With the exception of the Affortit (1970) data, relatively good agreement is observed in Figure 

2.2-2, which compares the values measured by Gibby (1974), Leibowitz (1972), and Affortit with 

predicted heat capacities for a MOX fuel, (U0.8, Pu0.2)O2. Since the Affortit results are known to 

be generally low in comparison to results from other investigators, the correlation is considered to 

be in good agreement with the experimental data. 
2.6 

200

300

400

500

600

700

800

300 800 1300 1800 2300 2800 3300

Temperature (K)

Specific Heat Capacity (J/kg*K)

Hein

Leibowitz

Gronvold

FRAP

Figure 2.2-1 Specific heat capacity of UO2 from three experimenters compared with the 

FRAPCON-3/FRAPTRAN/MATPRO (FRAP) correlation (solid line) for 

UO2. 
2.7 

200

250

300

350

400

450

500

550

600

650

700

300 800 1300 1800 2300 2800 3300

Temperature (K)

Specific Heat Capacity (J/kg*K)

Gibby

Leibowitz

Affortit

FRAP

Figure 2.2-2 Specific heat capacity of (U0.8Pu0.2)O2 from three experimenters compared 

with the FRAPCON-3/FRAPTRAN/MATPRO (FRAP) correlations (solid 

line) for MOX. 

2.2.2 Applicability and Uncertainty 

The correlations for fuel specific heat and enthalpy are valid for temperatures from 300K to more 

than 4000K. The standard error of the UO2 specific heat capacity correlation is  3 J/kg*K. 

However, this error increases at temperatures above 2300K, as seen in Figure 2.2-1, and may lead 

to underpredictions and overpredictions of 59 and 66 J/kg*K, respectively. As for the mixedoxide-specific heat capacity correlation, it is 6 to 10 J/kg*K, depending on the fraction of PuO2. 

For nonstoichiometric fuels, these uncertainties are approximately double. 

2.2.3 References 

Kerrisk, J.F., and D.G. Clifton. 1972. “Smoothed Values for the Enthalpy and Heat Capacity of 

UO2,” Nuclear Technology, 16:531-535. 

Kruger, O.L., and H. Savage. 1968. “Heat Capacity and Thermodynamic Properties of 

Plutonium Dioxide,” Journal of Chemical Physics, 45:4540-4544. 

Leibowitz, L., et al. 1971. “Enthalpy of Liquid Uranium Dioxide to 3500 K,” Journal of Nuclear 

Materials, 39:115-116. 

Hein, R.A., L.H. Sjodahl, and R. Swarc. 1968. “Heat Content of Uranium Dioxide from 1,200 to 

3,100 K,” Journal of Nuclear Materials, 25:99-102. 

 
2.8 

Leibowitz, L., L.W. Mschler, and M.G. Chasanov. 1969. “Enthalpy of Solid Uranium Dioxide 

from 2,500 K to its Melting Point,” Journal of Nuclear Materials, 29:356-358. 

Gronvold, F., et al. 1970. “Thermodynamics of the UO2+x Phase I. Heat Capacities of UO2.017

and UO2.254 from 300 to 1,000 K and Electronic Contributions,” Journal of Chemical 

Thermodynamics, 2:665-679. 

Affortit, C., and J. Marcon. 1970. “Chaleur Specifique a Haute Temperature des Oxydes 

D'Uranium et de Plutonium,” Revue Internationale des Hautes Temperatures et des Refractaires, 

7:236-241. 

Gibby, R.L., et al. 1974. “Analytical Expressions for Enthalpy and Heat Capacity for UraniumPlutonium Oxide,” Journal of Nuclear Materials, 50:155-161. 

Leibowitz, L., D.F. Fischer, and M.G. Chasanov. 1972. “Enthalpy of Uranium Plutonium 

Oxides (U0.8,Pu0.2) O1.07 from 2,350 to 3,000 K,” Journal of Nuclear Materials, 42:113-116. 

2.3 Fuel Thermal Conductivity (FTHCON) 

The subroutine FTHCON is used to calculate the thermal conductivity of the fuel pellet. Since 

accurate predictions of fuel rod behavior are strongly dependent on temperature, the calculation 

of fuel thermal conductivity is critical. The FTHCON subroutine is used in FRAPCON-3.4, 

FRAPTRAN 1.4, and MATPRO. However, the calculation of fuel thermal conductivity has 

evolved from the model originally proposed in MATPRO to the modified Nuclear Fuels 

Industries (NFI) and Duriez/modified NFI models used in the FRAPCON-3.4 and FRAPTRAN 

1.4 codes for UO2, UO2-Gd2O3, and MOX, respectively. An additional fuel thermal conductivity 

model for MOX developed at Halden is included as an option in FRAPCON-3.4. Each of the 

models is discussed in detail in the following sections along with the evolution of the MATPRO 

model and the applicability and uncertainty of the models. 

Following the initial description of the MATPRO fuel thermal conductivity model, an overview 

of subsequent thermal conductivity models is given. This historical summary of how thermal 

conductivity modeling has evolved provides insight to model development and previous 

motivations for change. For example, the MATPRO model was followed by a model by Lucuta 

that accounted for burnup degradation and was incorporated in FRAPCON-3. However, there 

were several issues with the Lucuta model, including its non-standard form and weak burnup 

degradation term, that led to its replacement by a modified model developed at NFI. Further 

refinements of this model were implemented by Pacific Northwest National Laboratory (PNNL) 

to provide a better fit to UO2 data. In addition, PNNL added refinements to include the Duriezderived modifications to the NFI model, which permit MOX fuels to be modeled. These updates 

were incorporated in FRAPCON-3.3 and continue to be used in FRAPCON-3.4. 

2.3.1 Model Development and Comparisons 

The original MATPRO model for fuel thermal conductivity is based on a mechanistic description 

of thermal conductivity, including lattice vibration (i.e., phonon) and electron-hole pair 

(i.e., electronic) contributions (Hagrman et al., 1981). The MATPRO thermal conductivity model 

is a function of temperature, density, O/M ratio, and plutonium content. Burnup is also included; 

however, it is only used to calculate the melt temperature. Unlike the models used in 
2.9 

FRAPCON-3.4 and FRAPTRAN 1.4, MATPRO does not account for the degradation of thermal 

conductivity with increasing burnup. 

The dominant contributor to thermal conductivity (i.e., phonon or electron) is dependent on 

temperature. Phonon contributions to thermal conductivity are dominant below 1500K, whereas 

electronic contributions are dominant above 2000K. In general, thermal conductivity decreases 

with increasing porosity, but this effect is significantly decreased above 1600K due to radiation 

effects and gas conductivity within the pores at these temperatures. At higher temperatures 

(> 3000K), estimates of liquid fuel thermal conductivity are made from physical considerations 

because data for molten fuel have not been found. 

The MATPRO model used to express thermal conductivity of solid fuel is described with 

Equation 2.3-1. 

(2.3-1) 

where 

k = thermal conductivity (W/m*K). 

D = fraction of theoretical density (unitless). 

CV = phonon contribution to the specific heat at constant volume (J/kg*K). The first 

term of the correlation for fuel specific heat capacity is used for this factor. 

eth = linear strain caused by thermal expansion when temperature is > 300K (unitless); 

the MATPRO correlation for fuel thermal expansion is used for this factor 

T = fuel temperature (K). 

T  = porosity correction for temperature, for temperature < 1364, T  = 6.50 – 

T*(4.69x10-3), for temperature > 1834, T  = -1, and for temperatures in the range 

from 1364 to 1834K, T  is found by interpolation. 

T  = fuel temperature if < 1800K. For temperatures > 2300K, T is equal to 2050K; 

for temperature in the range 1800 to 2300K, T is found by interpolation. 

A = a factor proportional to the point defect contribution to the phonon mean free 

path (m*s/kg*K); the correlation used for this factor is 0.339 + 12.6 x absolute 

value (2.0 – O/M ratio). 

B = a factor proportional to the phonon-phonon scattering contribution to the phonon 

mean free path (m*s/kg*K); the correlation used for this factor is 0.06867 x (1+ 

0.6238 x plutonium content of the fuel). 

The first term of Equation 2.3-1 represents the phonon contribution to specific heat. The 

temperature and density dependence of this term can be seen in Equation 2.3-1. Although not as 

evident, the dependence of the first term on O/M ratio and plutonium content is also included. 

The specific heat capacity of the fuel (CV) and the phonon scattering factor (B) are dependent on 

the plutonium content, while the point defect contribution factor (A) is dependent on the O/M 

ratio. Temperature dependence can be observed in the first and second terms of Equation 2.3-1, 

which represents the electronic contribution. 

     

          3 2

th

V

2T/13358169.01T/13358expT10x2997.5

e31''TBA

C

D1'T00469.05.61

D k

   

















 















  

 


2.10 

This MATPRO model (Hagrman et al., 1981) was initially replaced in FRAPCON-3 and 

FRAPTRAN with a model from Lucuta (Lucuta et al., 1996) which accounts for the degradation 

of thermal conductivity due to burnup. The Lucuta model for fuel thermal conductivity, K, is 

presented in Equation 2.3-2 

(2.3-2) 

where 

Ko = thermal conductivity of unirradiated, fully dense urania 

FD = factor for dissolved fission products 

FP = factor for precipitated fission products 

FM = factor to correct for the Maxwell porosity effect 

FR = factor for the radiation effect 

Thermal conductivity of unirradiated, fully dense urania and factors included in the Lucuta model 

are described by the Equations 2.3-3 through 2.3-7. 

(2.3-3) 

(2.3-4) 

(2.3-5) 

(2.3-6) 

(2.3-7) 

where 

Ko = conductivity of unirradiated, fully dense urania (W/m-K) 

T = temperature (K) 

B = burnup in atom% (1 atom% = 9.383 GWd/MTU at 200 MeV/fission) 

p = porosity fraction (as-fabricated plus swelling) 

s = shape factor ( = 1.5 for spherical pores) 

Although the Lucuta model for thermal conductivity accounted for burnup degradation, this 

model tended to overpredict thermal conductivity values at high temperature (> 2200K). 

 o     FRFMFPFDKK





 



  



 







   T

16361

exp T

10x715.4

T10x165.20375.0

1 K 2

9

o 4































 







T

B

0643.0

B

09.1

1 arctanT

B

0643.0

B

09.1 FD

265.3

265.3





















 





 





 

) 100

1200T exp(1

1

B019.03

B019.0 1FP

p)1s(1

p1 FM



 

) 80

900T exp(1

2.0 1FR  


2.11 

Furthermore, the model had too little burnup degradation, making it non-conservative above 

30 GWd/MTU. In addition, the form of the Lucuta model is non-standard, which does not 

facilitate comparison between models. As a result, the Lucuta model was later replaced with a 

modified version of a thermal conductivity model proposed by NFI (Ohira and Itagaki, 1997) 

(Lanning et al., 2005). 

 

Both FRAPCON-3.4 and FRAPTRAN 1.4 currently model urania fuel pellet thermal conductivity 

with the modified version of the pellet thermal conductivity model proposed by NFI (Ohira and 

Itagaki, 1997). The original NFI model was modified to alter the temperature-dependent portion 

of the burnup function in the phonon terms and change the electronic term (Lanning et al., 2005). 

The original (unmodified) and modified versions of the NFI thermal conductivity model are 

presented in Equations 2.3-8 and 2.3-9, respectively. 

(2.3-8) 

(2.3-9) 

where 

K = thermal conductivity, W/m-K 

T = temperature, K 

Bu = burnup, GWd/MTU 

f(Bu) = effect of fission products in crystal matrix (solution) 

 = 0.00187 * Bu 

g(Bu) = effect of irradiation defects 

 = 0.038 * Bu0.28

h(T) = temperature dependence of annealing on irradiation defects 

 = )T/Q(exp3961

1



Q = temperature-dependent parameter (“Q/R”) = 6380K 

A = 0.0452 m-K/W 

B = 2.46 x 10-4 m-K/W/K 

C = 5.47 x 10-9 W/m-K3

D = 2.29 x 1014 W/m-K5

E = 3.5 x 109

 W-K/m 

F = 16,361K 

The NFI model was further modified in FRAPCON-3.4 with a data correlation presented by 

Massih (Massih et al., 1992) to include gadolinia content. The modified NFI with gadolinia 

dependency used in FRAPCON-3.4 is presented in Equation 2.3-10. 

(2.3-10) 

2 4

95 DTCT

)T(h)Bu(g)Bu(fBTA

1 K 

 

)T/Fexp( T

E

)T(h)Bu(g))Bu04.0exp(9.01()Bu(fBTA

1 K95 2 

 

)T/Fexp( T

E

)T(h)Bu(g))Bu04.0exp(9.01()Bu(fBTgadaA

1 K

2 



 




2.12 

where 

a = constant = 1.1599 

gad = weight fraction of gadolinia 

As applied in FRACON-3.4 and FRAPTRAN 1.4, the above models are adjusted for as-fabricated 

fuel density (in fraction of theoretical density [TD]) using the Lucuta recommendation for 

spherical-shaped pores (Lucuta et al., 1996), as follows: 

(2.3-11) 

where 

Kd = thermal conductivity adjusted for as-fabricated fuel density, d 

K95 = thermal conductivity for 95 percent dense fuel 

d = as-fabricated fuel density 

The phonon-term modification in the modified NFI model reduces defect annealing at low burnup 

from that by Ohira and Itagaki (NFI). However, the original temperature-dependent annealing is 

restored at higher burnups such that, for burnups greater than 40 GWd/MTU, the phonon term is 

equivalent to that in the original NFI model. The electronic terms (which in either case become 

significant above 1500K) are altered in the modified NFI model to a more theoretically based 

equation. The magnitude is slightly lower than the original model at high temperature. This 

adjustment was indicated by the Institute for Transuranium Elements data on unirradiated 

pressurized-water reactor (PWR) pellet material at temperatures approaching fuel melting 

(Ronchi et al., 1999). 

At low burnups (< 20 GWd/MTU) and low temperatures (< 1000K), the modified model is higher 

than the unmodified NFI model and roughly equivalent to the Lucuta model without its radiation 

term. At higher burnup (> 30 GWd/MTU), the modified model is equivalent to the original NFI 

model with the exception of the small reduction at very high temperatures. 

Fuel thermal conductivity modeling was extended to include MOX in FRAPCON-3 and 

FRAPTRAN by adopting a model that was a combination of the Duriez stoichiometry-dependent 

correlation, which is derived from diffusivity measurements on unirradiated fuel pellets, and the 

modified version of the NFI model, which includes burnup degradation effects. The Duriez data 

places greater dependence of MOX thermal conductivity on stoichiometry, but only a minor 

dependence on plutonia content (Duriez et al., 2000). The Duriez/modified NFI model is 

presented in Equation 2.3-12. 

(2.3-12) 

where 

K95 = conductivity at 95 percent TD, W/m-K 

x = 2.00 – O/M (i.e., oxygen-to-metal ratio) 

T = temperature, K 

A(x) = 2.85x + 0.035, m-K/W 

d  95   )}]d1(5.00.1/{d[*K*0789.1K

) T

D exp( TC

)T(h)Bu(g))Bu4.0exp(9.01()Bu(fT)x(Bgada)x(A

1 K

2

mod

95

 

 




2.13 

B(x) = (2.86 – 7.15x)*1E-4 m/W 

C = 1.689x109

, W-K/m 

D = 13,520K 

a = 1.1599 

gad = weight fraction gadolinia (not expected in MOX) 

Bu = burnup in GWd/tHM 

f(Bu) = effect of fission products in crystal matrix (solution) 

 = 0.00187*Bu 

g(Bu) = effect of irradiation defects 

 = 0.038*Bu0.28

h(T) = temperature dependence of annealing on irradiation defects 

 = )T/Q(exp3961

1



Q = temperature dependent parameter (“Q/R”) = 6380K 

Cmod = 1.5x109

 W-K/m 

FRAPCON 3 provides an option to use a fuel thermal conductivity model derived at Halden. The 

Halden correlation does not include dependency on stoichiometry, but does emphasize burnup 

effects. The Halden correlation tends to overpredict thermal conductivity data at high 

temperatures (> 2500K) (Lanning et al., 2005). 

(2.3-13) 

where 

K95 = conductivity at 95 percent TD, W/m-K 

T = temperature, °C 

a = 1.1599 

gad = weight fraction gadolinia (not expected in MOX) 

B = burnup in MWd/kg UO2

 = minimum of 1650°C or current temperature in °C 

To summarize, the modified NFI model is used for UO2, whereas the Duriez/NFI model is used 

for MOX. The modified NFI model used in FRAPCON-3.4 accounts for gad, but the modified 

NFI model in FRAPTRAN 1.4 does not. The Duriez/NFI models account for gad in both codes, 

although no gad additions are included in MOX at this time. Finally, FRAPCON-3.4 has the 

option for Halden thermal conductivity correlation, but FRAPTRAN does not. 

A comparison between the MATPRO model and the modified NFI model is presented in Figure 

2.3.3-1. The MATPRO model predicts lower thermal conductivities than the modified NFI 

model across the applicable temperature range (i.e., 500 to 3000K) for unirradiated UO2. 

However, unlike the MATPRO model, the modified NFI model can account for degradation of 

thermal conductivity with increasing burnup. The modified NFI model at a burnup level of 

30 GWd/MTU is shown for comparison and reveals significant degradation of thermal 

conductivity. This represents an improvement relative to the MATPRO model, which does not 

account for burnup degradation of thermal conductivity. 

)T00188.0exp(0132.0

)B00333.01(10x475.2B0040.0x1599.1gada1148.0

92.0 K95 4



    




2.14 

0.00

1.00

2.00

3.00

4.00

5.00

6.00

7.00

8.00

9.00

10.00

0 500 1000 1500 2000 2500 3000

Temperature (K)

Thermal Conductivity (W/m*K)

MATPRO

Modified NFI -

0 GWd/MTU

Modified NFI -

30GWdMTU

Figure 2.3.3-1 MATPRO and modified NFI thermal conductivity models for UO2 as a 

function of temperature. Unlike the MATPRO model, the modified NFI 

model can account for burnup and is presented at burnup levels of 0 and 

30 GWd/MTU for comparison. 

Model-to-data comparisons for the MATPRO and modified NFI models and data collected from 

unirradiated UO2 are made in Figures 2.3.3-2 through 2.3.3-3. In addition to Ronchi (1999), the 

unirradiated UO2 data sets from that were used to derive the MATPRO model are also included. 

Based on Figure 2.3.3-2, it can be seen that the MATPRO model tends to underpredict the 

thermal conductivity of the unirradiated UO2. The average bias in Figure 2.3.3-2 underpredicts 

thermal conductivity by about 9.5 percent, with a standard deviation of about 11 percent. The 

modified NFI model also underpredicts the unirradiated data, but only by about 2.4 percent and 

with a standard deviation of 8 percent (see Figure 2.3.3-3). Based on the smaller bias and 

standard deviation, it can be seen that the modified NFI model predicts the thermal conductivity 

of UO2 much better than the MATPRO model. 
2.15 

MATPRO (UO2)

1

2

3

4

5

6

7

8

12345678

Measured Thermal Conductivity (W/m*K)

Predicted Thermal 

Conductivity (W/m*K)

Predicted = Measured Christensen Godfrey

Bates Gibby Weilbacher

Goldsmith and Douglas Hobson Ronchi

Figure 2.3.3-2 MATPRO model-to-data comparison for thermal conductivity of unirradiated 

UO2. 

Modified NFI (UO2)

1

2

3

4

5

6

7

8

12345678

Measured Thermal Conductivity (W/m*K)

Predicted Thermal 

Conductivity (W/m*K)

Predicted = Measured Christensen Godfrey

Bates Gibby Weilbacher

Goldsmith and Douglas Hobson Ronchi

Figure 2.3.3-3 Modified NFI model-to-data comparison for thermal conductivity of 

unirradiated UO2. 
2.16 

Model-to-data comparisons for the MATPRO and modified NFI models and data collected from 

irradiated UO2 are made in Figures 2.3.3-4 through 2.3.3-5. Based on Figure 2.3.3-4, it can be 

seen that the MATPRO model predicts significantly higher values for thermal conductivity than 

the measured values obtained from irradiated UO2. The average bias in Figure 2.3.3-4 indicates 

that the thermal conductivity is overpredicted by about 40 percent, with a standard deviation of 

about 31 percent. The modified NFI model also tends to overpredict the small amount of 

irradiated data on average, but only by about 5 percent and with a standard deviation of about 

8 percent (see Figure 2.3.3-5). Based on the significantly smaller bias and standard deviation, it 

can be seen that the modified NFI model predicts the thermal conductivity of irradiated UO2

better than the MATPRO model. The disparity between models is attributed to the degradation of 

thermal conductivity with increasing burnup, which is accounted for in the modified NFI model 

but not the MATPRO model. 

1

2

3

4

5

12345

Predicted Thermal Conductivity (W/m*K)

Measured Thermal Conductivity (W/m*K)

MATPRO (UO2)

Predicted = Measured Ronchi 30-35 GWd/MTU Carrol 40 GWd/MTU

Ronchi 50-55 GWd/MTU Ronchi 82 GWd/MTU Ronchi 92 GWd/MTU

Ronchi 96 GWd/MTU

Figure 2.3.3-4 MATPRO model-to-data comparison for thermal conductivity of irradiated 

UO2. 
2.17 

Modified NFI (UO2)

1

2

3

4

5

12345

Measured Thermal Conductivity (W/m*K)

Predicted Thermal 

Conductivity (W/m*K)

Predicted = Measured Ronch 30-35 GWd/MTU Carrol 40 GWd/MTU

Ronchi 50-55 GWd/MTU Ronchi 82 GWd/MTU Ronchi 92 GWd/MTU

Ronchi 96 GWd/MTU

Figure 2.3.3-5 Modified NFI model-to-data comparison for thermal conductivity of 

irradiated UO2. 

A comparison between the MATPRO model and the Duriez/modified NFI model is presented in 

Figure 2.3.3-6 for MOX fuel with a plutonium concentration of 7 wt%. The MATPRO model 

predicts lower thermal conductivities than the Duriez/modified NFI model across the applicable 

temperature range (i.e., 500 to 3000K). However, unlike the MATPRO model, the 

Duriez/modified NFI model can account for degradation of thermal conductivity with increasing 

burnup. The Duriez/modified NFI model at a burnup level of 30 GWd/MTU is shown for 

comparison and reveals significant degradation of thermal conductivity. This represents an 

improvement relative to the MATPRO model, which does not account for burnup degradation of 

thermal conductivity. 
2.18 

0.00

1.00

2.00

3.00

4.00

5.00

6.00

7.00

8.00

9.00

10.00

0 500 1000 1500 2000 2500 3000

Temperature (K)

Thermal Conductivity (W/m*K)

MATPRO-MOX

Duriez-Modified NFI

Duriez-Modified NFI

30GWd/MTU

Figure 2.3.3-6 MATPRO and Duriez/modified NFI thermal conductivity models as a 

function of temperature for MOX fuel with 7 wt% plutonium content. 

Unlike the MATPRO model, the Duriez/modified NFI model can account for 

burnup and is presented at burnup levels of 0 and 30 GWd/MTU for 

comparison. 

Model-to-data comparisons for the MATPRO and Duriez/modified NFI models and unirradiated 

MOX fuel are made in Figures 2.3.3-7 and 2.3.3-8. The comparisons utilize the data used to 

derive the Duriez/modified NFI model. Data collected from MOX with plutonium concentrations 

of 3 and 15 wt% were used (Duriez et al., 2000). Based on Figure 2.3.3-7, it can be seen that the 

MATPRO model tends to underpredict thermal conductivity. The average bias in Figure 2.3.3-7 

underpredicts thermal conductivity by about 9 percent at both 3 and 15 wt% plutonium 

concentrations. The standard deviation of the MATPRO predictive bias at both plutonium 

concentrations is under 6 percent. As expected, the Duriez/modified NFI model exhibits better 

agreement with the limited experimental data from which it was derived. The average bias 

indicates that the MOX model overpredicts by less than 2 percent and the standard deviation 

between the model and the data is less than 3 percent. Based on the smaller bias and standard 

deviation, it can be seen that the Duriez/modified NFI model fits this MOX thermal conductivity 

data better than the MATPRO model. 

It should be noted that the MOX fuel thermal conductivity is strongly influenced by O/M ratio, 

with conductivity decreasing with a deviation from stoichiometry (2.0). Based on Figure 2.3.3-8, 

it can be seen that plutonium content has a relatively minor effect at concentrations less than 

15 wt%. The current MOX fuel is fabricated closer to being stoichiometric than fuel typically 

fabricated more than 25 years ago. The older data from which the MATPRO model was 

developed were not as close to being stoichiometric as the more recent data. 
2.19 

MATPRO (MOX)

1

2

3

4

5

12345

Measured Thermal Conductivity (W/m*K)

Predicted Thermal 

Conductivity (W/m*K)

Predicted = Measured 3 wt% Pu 15 wt% Pu

Figure 2.3.3-7 MATPRO model-to-data comparison for thermal conductivity of MOX 

modeled with 7 wt% of plutonium. 

Duriez-Modified NFI (MOX)

1

2

3

4

5

12345

Measured Thermal Conductivity (W/m*K)

Predicted Thermal 

Conductivity (W/m*K)

Predicted = Measured 3 wt% Pu 15 wt% Pu

Figure 2.3.3-8 Duriez/modified NFI model-to-data comparison for thermal conductivity of 

MOX with modeled with 7 wt% of plutonium. 
2.20 

2.3.2 Applicability and Uncertainty 

The following are limits that restrict the applicability of the MATPRO fuel thermal conductivity 

model (Equation 2.3-1): 

 Theoretical density must range between 90 and 100 percent. 

 Temperature must range between 500 (the Debye temperature) and 3000K. 

 Above 3000K, it is assumed that the fuel becomes molten and the uncertainties on predictions 

made above this temperature are large. 

In addition, several assumptions have been made to provide an approximate model for effects of 

variation in the plutonium content and O/M ratio of ceramic fuels: 

 The effect of variation in theoretical density for MOX fuels has been assumed to be described 

by the porosity correction derived with UO2 data. 

 The high-temperature electronic contribution to thermal conductivity has been assumed to be 

the same for PuO2, UO2, UO2-Gd2O3, and nonstoichiometric fuels. 

 Variations in plutonium content have been assumed to affect only the phonon-phonon 

scattering factor. 

 Variation in O/M ratio has been assumed to affect only the defect term. 

The bias and relative error of the MATPRO fuel thermal conductivity model was assessed based 

on the model-to-data comparisons previously presented for unirradiated UO2, irradiated UO2, and 

MOX. The average and standard deviations of the bias for unirradiated UO2 are approximately 

-9.5 and 11 percent, respectively, while the relative error is approximately 19 percent. For 

irradiated UO2, which has a limited number of data points for comparison, the average and 

standard deviations of the bias are 40 and 31 percent, respectively. The relative error for the 

limited number of irradiated UO2 data points is approximately 51 percent. Finally, for MOX, the 

average and standard deviations of the bias are -9 and 6 percent, respectively, while the relative 

error is approximately 13 percent. For reference, a negative bias indicates that the model 

underpredicts the data, whereas a positive bias indicates that the model overpredicts the data 

For the modified NFI and the Duriez/modified NFI models, the recommended ranges for 

application for modeling the thermal conductivity of UO2 are as follows: 

 Temperature: 300 to 3000K 

 Rod-average burnup: 0 to 62 GWd/MTU 

 As-fabricated density: 92 to 97 percent TD 

 Gadolinia content: 0 to 10 wt% 

For the Duriez/modified NFI model, the thermal conductivity model is applicable for plutonia 

particle sizes of less than 20 microns. 

The bias and relative error of the modified NFI and the Duriez/modified NFI fuel thermal 

conductivity models were assessed based on the model-to-data comparisons previously presented 

for unirradiated UO2, irradiated UO2, and MOX. The average and standard deviations of the bias 

for unirradiated UO2 are approximately -2.4 and 8 percent, respectively, while the relative error is 

approximately 9 percent. For irradiated UO2, which has a limited number of data points for 
2.21 

comparison, the average and standard deviations of the bias are approximately 5 and 8 percent, 

respectively, while the relative error is approximately 10 percent. Finally, for the 

Duriez/modified NFI fuel thermal conductivity model for MOX, the average and standard 

deviations of the bias are less than 2 and 3 percent, respectively, while the relative error is 

approximately 3 percent. These are values are based on comparisons to the small amount of nearstoichiometric MOX data. For reference, a negative bias indicates that the model underpredicts 

the data, whereas a positive bias indicates that the model overpredicts the data 

2.3.3 References 

Hagrman, D.L., G.A. Reymann, and G.E. Mason. 1981. A Handbook of Materials Properties for 

Use in the Analysis of Light Water Reactor Fuel Rod Behavior, MATPRO Version 11 (Revision 

2), NUREG/CR-0479 (TREE-1280), prepared by EG&G Idaho, Inc., Idaho Falls, ID, for the U.S. 

Nuclear Regulatory Commission, Washington, D.C. 

Lucuta, P.G., H.S. Matzke, and I.J. Hastings. 1996. “A Pragmatic Approach to Modeling 

Thermal Conductivity of Irradiated UO2 Fuel: Review and Recommendations,” Journal of 

Nuclear Materials, 232:166-180. 

Ohira, K., and N. Itagaki. 1997. “Thermal Conductivity Measurements of High Burnup UO2

Pellet and a Benchmark Calculation of Fuel Center Temperature,” in Proceedings of the ANS 

International Topical Meeting on LWR Fuel Performance, Portland, Oregon, March 2-6, 1997, 

pp. 541-549. 

Lanning, D.D., C.E. Beyer, and K.G. Geelhood. 2005. FRAPCON-3 Updates, Including MixedOxide Fuel Properties, NUREG/CR-6534, Vol. 4, PNNL-11513, Pacific Northwest National 

Laboratory, Richland, WA. 

Massih, A.R., S. Persson, and Z. Weiss. 1992. “Modeling of (U,Gd)O2 Fuel Behavior in Boiling 

Water Reactors,” in Proceedings of Symposium E on Nuclear Materials for Fission Reactors of 

the 1991 E-MRS Fall Conference; Journal of Nuclear Materials, 188:323-330. 

Ronchi, C., M. Sheindlin, M. Musella, and G.J. Hyland. 1999. “Thermal Conductivity of 

Uranium Dioxide Up to 2900K from Simultaneous Measurement of the Heat Capacity and 

Thermal Diffusivity,” Journal of Applied Physics, 85(2):776-789. 

Duriez, C., J.-P. Allesandri, T. Gervais, and Y. Philipponnea. 2000. “Thermal Conductivity of 

Hypostoichiometric Low Pu Content (U,Pu)O2-x Mixed Oxide,” Journal of Nuclear Materials, 

277:143-158. 

Bates, J.L., C.A. Hinman, and T. Kawada. 1967. “Electrical Conductivity of Uranium Dioxide,” 

Journal of the American Ceramic Society, 50:652-656. 

Goldsmith, L.A., and J.A.M. Douglas. 1973. “Measurements of the Thermal Conductivity of 

Uranium Dioxide at 670-1270 K,” Journal of Nuclear Materials, 47:31-42. 

Christensen, J.A. et al. 1964. “Uranium Dioxide Thermal Conductivity,” Transactions of the 

American Nuclear Society, 7:391-392. 
2.22 

Gibby, R.L. 1971. “The Effect of Plutonium Content on the Thermal Conductivity of (U,Pu)O2

Solid Solutions,” Journal of Nuclear Materials, 38:163-177. 

Hobson, I.C., R. Taylor, and J.B. Ainscough. 1974. “Effect of Porosity and Stoichiometry on the 

Thermal Conductivity of Uranium Dioxide,” Journal of Physics, Section D: Applied Physics, 

7:1003-1015. 

Godfrey, T.G., et al. 1964. Thermal Conductivity of Uranium Dioxide and Armco Iron by an 

Improved Radial Heat Flow Technique, ORNL-3556, Oak Ridge National Laboratory, Oak 

Ridge, TN. 

Weilbacher, J.C. 1972. “Diffusivite Thermique de l’Oxyde d’Uranium et de l’Oxyde de 

Thorium a Haute Temperature,” High Temperatures-High Pressure, 4:431-438. 

Ronchi, C., M. Sheindlin, D. Staicu, and M. Kinoshita. 2004. “Effect of Burn-up on the Thermal 

Conductivity of Uranium Dioxide up to 100,000 MWdt-1,” Journal of Nuclear Materials, 327:58-

76. 

2.4 Fuel Emissivity (FEMISS) 

The subroutine FEMISS is used to calculate the total hemisphericalfuel emissivity (emissivity 

integrated over all wavelengths) as a function of temperature. Fuel emissivity is defined as the 

ratio of radiant energy emitted from a material to that emitted by a black body at the same 

temperature. The subroutine FEMISS is used to calculate radiant energy transfer from fuel to 

cladding in conjunction with thermal conduction. Radiant energy transfer can be a significant 

heat transfer mechanism, depending on the gap size, temperature gradient across the gap, and 

plenum gas. The FEMISS subroutine used by FRAPCON-3 and FRAPTRAN is the same as the 

subroutine documented in MATPRO. 

2.4.1 Model Development and Comparisons 

According to the Stefan-Boltzmann law, the total radiant power per unit area emitted by a body at 

temperature T is 

(2.4-1) 

where 

P = radiant power per unit area (W/m2

) 

e = total hemispherical emissivity (unitless) 

 = the Stefan-Boltzmann constant (5.672 x 10-8 W/m2

-K) 

T = temperature (K) 

The expression used in the FEMISS subroutine to describe total emissivity is 

(2.4-2) 

The first term of Equation 2.4-2 is slightly less than the value used in MATPRO, 0.7856. The 

standard error of Equation 2.4-2 with respect to its database is  6.796 percent, which is in 

4  TeP

T10x5263.178557.0e 5 
2.23 

agreement with the error reported in MATPRO,  6.8 percent. These small differences are 

attributed to rounding error and do not significantly alter emissivity predictions. 

The emissivity data were measured at temperatures up to approximately 2400K, and use of 

FEMISS above this temperature is speculative because of possible high-temperature effects that 

are not modeled. At the time of model development, there were no data to develop a MOX 

emissivity equation, so Equation 2.4-2 is also recommended for MOX. 

The hemispherical spectral data of Held and Wilder (1969) and the emissivity data of Cabannes 

(1967) and J.M. Jones and D.G. Murchison (1965) were used in developing the FEMISS model. 

The emissivity data used to model comparison for FEMISS are presented in Figure 2.4-1. 

0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

300 800 1300 1800 2300 2800

Temperature (K)

Emissivity (unitless)

Held and Wilder

Cabannes et al.

FEMISS

Standard Error

Figure 2.4-1 Data to model comparison for FEMISS. 

2.4.2 Applicability and Uncertainty 

The correlation used in the FEMISS subroutine to calculate fuel emissivity is applicable at 

temperatures up to 2400K. The standard error of this correlation with respect to the database 

from which it was derived is approximately  6.8 percent 

2.4.3 References 

Held, P.C., and D.R. Wilder. 1969. “High Temperature Hemispherical Spectral Emittance of 

Uranium Oxides at 0.65 and 0.70 m,” Journal of the American Ceramic Society, 52(4):182-185. 

M.M.F. Cabannes, J.P. Stora, and J. Tsakiris. 1967. “ORTIQUE-MOLECULAIRE-Fracteurs de 

re’flexion et d’e’mission de UO2 a’ haute Temperature,” C.R. Acad. SC. Paris, 264. 

J.M. Jones and D.G. Murchison. 1965. “Optical Properties of Uranium Oxides,” Nature, 

205(4972):663-665. 
2.24 

2.5 Fuel Thermal Expansion (FTHEXP) 

The subroutine FTHEXP models dimensional changes in unirradiated fuel pellets caused by 

thermal expansion. It can handle any combination of UO2, UO2-Gd2O3, or PuO2 in solid, liquid, 

or solid-liquid states and includes expansion due to the solid-liquid phase change. Dimensional 

changes in the fuel affect the pellet-to-cladding gap size, which is a major factor in determining 

gap heat transfer and thus the stored energy, an important quantity for safety analysis. The 

subroutine FTHEXP is used in FRAPCON-3.4, FRAPTRAN 1.4, and MATPRO. However, the 

FRAPCON-3 and FRAPTRAN codes use updated fitting coefficients based on more recent data. 

2.5.1 Model Development and Comparisons 

The FTHEXP subcode models fuel thermal expansion as a function of temperature, fraction of 

PuO2, and the fraction of fuel which is molten. The O/M ratio is not included. When the 

departure from stoichiometry (O/M – 2.0) is greater than 0.2, there is clearly an effect. However, 

this effect is ignored in modeling thermal expansion, since typical reactor fuels only deviate less 

than a tenth this much from the stoichiometric composition. 

The equations for thermal expansion of UO2, UO2-Gd2O3, and PuO2 have the same form. In the 

solid phase, Equation 2.5-1 is used. 

 (2.5-1) 

where 

L/L = linear strain caused by thermal expansion (equal to zero at 300K) (unitless) 

T = temperature (K) 

ED = energy of formation of a defect (J) 

k = Boltzmann’s constant (1.38 x 10-23 J/K) 

and K1, K2, and K3 are constants to be determined. 

The correlation given in Equation 2.5-1 was originally used in MATPRO and was incorporated 

into both the FRAPCON-3 and FRAPTRAN codes. Data collected by Baldock et al. (1966), 

Burdick and Paker (1956), Gronvold (1955), Christensen (1963), and Kempeter and Elliott (1958) 

were used to determine the correlation constants for UO2 used in MATPRO. However, newer 

data provided by Martin (1988) and Momin et al. (1991) required the constants to be updated to 

improve the fit between the correlation and high-temperature data. These updated constants are 

included in FRAPCON-3.4/FRAPTRAN 1.4 (FRAP). Table 2.5-1 presents the original 

MATPRO constants and the updated constants for UO2, as well as the constants for PuO2. 

Table 2.5-1 Parameters used in UO2 and PuO2 solid-phase thermal expansion 

correlations. 

Constant 

MATPRO 

UO2

FRAP 

UO2 PuO2 Units 

K1 1.0 x 10-5 9.80 x 10-6 9.0 x 10-6 K-1 

K2 3.0 x 10-3 2.61 x 10-3 2.7 x 10-3 Unitless 

K3 4.0 x 10-2 3.16 x 10-1 7.0 x 10-2 Unitless 

ED 6.9 x 10-20 1.32 x 10-19 7.0 x 10-20 J 

  10   32  D )kT/Eexp(KKTKL/L
2.25 

For mixed UO2 and PuO2, the thermal expansion of the solid is found by combining the 

contribution from each constituent in proportion to its weight fraction. 

During melting, an expansion equal to a linear strain of 0.043 occurs. If the fuel is partially 

molten, the strain due to thermal expansion is given by Equation 2.5-2. 

 (2.5-2) 

where 

L/L0(Tm) = thermal expansion strain of solid fuel from equations with T = Tm

Tm = melting temperature of the fuel (K) 

FACMOT = fraction of the fuel which is molten (unitless) 

 If FACMOT = 0.0, the fuel is all solid; 

 If FACMOT = 1.0, the fuel is all molten 

The correlation used to describe the expansion of entirely molten fuel is given by Equation 2.5-3. 

(2.5-3) 

The solid-to-liquid phase transition is isothermal only for pure UO2 or pure PuO2. For MOX, the 

transition occurs over a finite temperature range, denoted in Equation 2.5-3 by Tm. 

Comparisons between the databases for UO2 and PuO2 with their respective correlations are 

presented in Figures 2.5-1 and 2.5-2, respectively. In Figure 2.5-1, it can be seen that the updated 

constants used in FRAPCON-3.4/FRAPTRAN 1.4 (FRAP) for UO2 provide a better fit to the 

high-temperature data than the previous constants used in MATPRO. The comparison in Figure 

2.5-2 reveals good agreement between the correlation and the database. 

0.000

0.005

0.010

0.015

0.020

0.025

0.030

0.035

0.040

0.045

0.050

0 500 1000 1500 2000 2500 3000

Temperature (K)

Fractional Length Change (L/L)

Baldock Burdick Gronvold

Other MATPRO References Martin References MATPRO UO2

FRAP UO2

Figure 2.5-1 Model-to-data comparison for UO2 correlation used in MATPRO and 

FRAPCON-3/FRAPTRAN with updated constants. 

0   m0  043.0)T(L/LL/L  FACMOT

   m m

5

0  m0   TTT10x6.3043.0)T(L/LL/L 
2.26 

0.000

0.002

0.004

0.006

0.008

0.010

0.012

0.014

0.016

0 500 1000 1500 2000

Temperature (K)

Fractional Length Change (L/L)

FRAP and MATPRO PuO2 Tokar and Brett References

Figure 2.5-2 Model-to-data comparison for PuO2 correlation used in MATPRO and 

FRAPCON-3/FRAPTRAN. 

2.5.2 Applicability and Uncertainty 

The fuel thermal expansion correlation used in the FTHEXP subroutine is applicable through the 

melting temperature of the fuel (3111.15K for UO2). The uncertainty of the pooled data was 

found to be temperature dependent, increasing approximately linearly with temperature. 

Therefore, a percentage error is given rather than a fixed number. The   limits were found to 

be within 10 percent of the calculated value. 

2.5.3 References 

Baldock, P.J., et al. 1966. “The X-ray Thermal Expansion of Near-Stoichiometric UO2,” Journal 

of Nuclear Materials, 18:305-313. 

Burdick, M.D., and H.S. Parker. 1956. “Effect of Particle Size on Bulk Density and Strength 

Properties of Uranium Dioxide Specimens,” Journal of the American Ceramic Society, 39:181-

187. 

Gronvold, F. 1955. “High Temperature X-ray Study of Uranium Oxides in the UO2-U3O8

Region,” Journal of Inorganic and Nuclear Chemistry, 1:357-370. 

J.A. Christensen. 1963. “Thermal Expansion and Change in Volume of Uranium Dioxide on 

Melting” Journal of the American Ceramic Society, 46:607-608. 

Kempter, C.P., and R.E. Elliott. 1958. “Thermal Expansion of UN, UO2, UO2-ThO2, and ThO2,” 

Journal of Chemical Physics, 30:1524-1526. 
2.27 

Martin, D.G. 1988. “The Thermal Expansion of Solid UO2 and (U,Pu) Mixed Oxides – A 

Review and Recommendations,” Journal of Nuclear Materials, 152:94-101. 

Momin, A.C., E.B. Mirza, and M.D. Mathews. 1991. “High Temperature X-Ray 

Diffacractometric Studies on the Lattice Thermal Expansion Behavior of UO2, ThO2, and 

(U0.2Th0.8)O2 Doped with Fission Products,” Journal of Nuclear Materials, 185:308-310. 

Tokar, M., and A.W. Nutt. 1972. “Thermal Expansion of PuO2 from 25 to 1420°C,” 

Transactions of the American Nuclear Society, 10:210-211. 

Brett, N.H., and L.E. Russel. 1960. “The Thermal Expansion of PuO2 and Some Other Actinide 

Oxides Between Room Temperature an 1000°C,” in Proceedings of the Second International 

Conference on Plutonium Metallurgy, Grenoble, France, April 19-22, 1960, pp. 397-410. 

2.6 Fuel Densification (FUDENS) 

The subroutine FUDENS calculates fuel dimensional changes due to densification of UO2 and 

MOX during the first few thousand hours of water reactor operation. If available, data from 

resintering tests should be input by the user to determine densification (Regulatory Guide 1.126 

[NRC 1978]). Otherwise, the subcode FUDENS calculates densification as a function of fuel 

burnup, temperature, and initial density. 

The data used to develop FUDENS were taken from irradiated fuel which also included swelling 

(Rolstad, 1974). If fuel densification is much greater than swelling during the first 1000 hours of 

irradiation, then, to a first approximation, swelling can be neglected during this period. This was 

done during the development of the FUDENS model. 

The FUDENS subroutine is used in FRAPCON-3.4 and is similar to the correlation described by 

MATPRO. However, FUDENS is not included in FRAPTRAN because FRAPTRAN is intended 

for transient events occurring over a short time scale. Since densification occurs over longer time 

scales, it is not included in FRAPTRAN. 

2.6.1 Model Development and Comparisons 

The subroutine FUDENS uses one of two methods to calculate the maximum density change 

during irradiation. The RSNTR method uses the density change observed during resintering tests 

(1973K for 24 hours based on Regulatory Guide 1.126 [NRC 1978]) in a laboratory furnace and 

is the preferred input for the calculation. When RSNTR > 0, MATPRO uses either 

Equation 2.6-1 or 2.6-2 and 2.6-6 while FRAPCON-3.4 uses Equations 2.6-3 and 2.6-6. If a 

resintering density change is not input (i.e., RSNTR = 0), the TSINT method uses the initial 

unirradiated density of the fuel and the fuel fabrication sintering temperature and burnup for 

density calculations. For TSNTR calculations, both MATPRO and FRAPCON-3.4 use either 

Equation 2.6-4 or 2.6-5 and 2.6-6. Both the RSNTR and TSNTR methods for MATPRO and 

FRAPCON-3.4 are described by the following equations to calculate the maximum densification 

length change during irradiation. 

In the MATPRO version of FUDENS, if a nonzero value for the resintering density change is 

input, then either Equation 2.6-1 or Equation 2.6-2 is used, depending on the fuel temperature, 

FTEMP (Siefken et al., 2001). 
2.28 

For FTEMP < 1000K (2.6-1) 

For FTEMP  1000K (2.6-2) 

In FRAPCON-3.4, if the resintering density is a nonzero input, then the following equation is 

used to calculate the maximum dimensional change. 

(2.6-3) 

If zero is input for the resintering density change, then the MATPRO and FRAPCON-3.4 models 

use either Equation 2.6-4 or Equation 2.6-5, depending on the fuel temperature, FTEMP (Rolstad 

et al., 1974). 

For FTEMP < 1000K (2.6-4) 

For FTEMP  1000K (2.6-5) 

where 

L m

L 



 



 

 = maximum possible dimension change of fuel due to irradiation (percent) 

RSNTR = resintered fuel density change (kg/m3

) 

FTEMP = fuel temperature (K) 

DENS = initial density (percent theoretical) 

TSINT = sintering temperature (K) 

FDENS = starting density (kg/m3

) 

Densification as a function of burnup is calculated using 

(2.6-6) 

where 





 



 

L

L

 = dimension change (percent) 

FBU = fuel burnup (MWd/kgU) 

  0015.0 RSNTR

L

L

m

  

 



 

  00285.0 RSNTR

L

L

m

  

 



 

 

  TSINT 1453

1002.22 DENS

L

L

m 

    

 



 

 

  TSINT 1453

1006.66 DENS

L

L

m 

    

 



 

          BFBU3 BFBU35

m

e e0.2

L

L

L

L      

 



   

)FDENS*0.3/(RSNTR*100 L

L

m

  

 



 
2.29 

B = a constant determined by the subcode to fit the boundary condition: L/L 

= 0 when FBU = 0 

The FUDENS subcode presented in Equation 2.6-6 is used to calculate total densification, and 

then the densification from the previous time step is subtracted to obtain the incremental 

densification. The incremental densification for the time step being considered is the output of 

the subcode FUDENS and is determined by Newton’s method (Hamming, 1971). 

The relationship between densification and burnup in Equation 2.6-6 was originally suggested by 

Rolstad et al. (1974). This relationship was adopted for use in the FUDENS subroutine because it 

successfully describes the burnup dependence of both the original data (Rolstad et al., 1974) and 

EPRI data (Freshley et al., 1976) that was collected shortly thereafter. Since these initial studies, 

there has been additional in-reactor densification data collected for UO2 and MOX fuels 

irradiated in thermal reactors. 

Banks (1974) reported densification data on solid and annular UO2 pellets irradiated in the 

Winfrith Steam Generating Heavy Water Reactor and the Halden Boiling Water Reactor. The 

solid pellets experienced centerline temperatures of 733 to 2443K over the burnup range 0.37 to 

23.7 MWd/kgU. They concluded that pellets with centerline temperatures less than 1423K 

densified during irradiation, but at higher temperatures dimensional change was affected by 

intergranular fission gas swelling. Therefore, only data for temperatures less than 1423K were 

included in comparison with the FUDENS model (Equations 2.6-4, 2.6-5, and 2.6-6). The 

mechanism for densification was determined to be the loss of pores with initial diameters less 

than 2 m. 

Freshley et al. (1978) reported densification data on MOX pellets irradiated in the radially 

adjustable facility tubes of the General Electric Test Reactor. Two pins were irradiated with fuel 

centerline temperatures up to 2048K and burnup to 4.6 MWd/kgM. The pellets were fabricated 

with sintering temperatures from 1748 to 2008K and had starting density values ranging from 91 

to 95 percent TD. They concluded that dilute additions of PuO2 in UO2 (up to 6 wt% PuO2) do 

not change the densification behavior of the fuel. The mechanism of densification was observed 

to be sintering of pores with initial diameters of 0.5 m or smaller. 

Small (1987) irradiated four types of UO2 pellets in the Harwell Materials Test Reactor. Pellet 

centerline temperatures ranged from 623 to 1823K and burnup ranged from 0.0025 to 

2.6 MWd/kgU. Three of the four types of pellets were fabricated with pores smaller than 2 m in 

diameter to more carefully examine the mechanisms behind densification. For reference, fuel 

currently fabricated has very few pores less than 2 m in diameter. Densification was determined 

by measuring porosity before and after irradiation. At centerline temperatures above 1273K, 

densification was caused by sintering of pores less than 1 m in diameter. 

Figure 2.6-1 compares measured and predicted values of densification using the TSNTR 

correlation, which is used in both FRAPCON-3.4 and MATPRO and incorporates Equations 

2.6-4, 2.6-5, and 2.6-6. 
2.30 

-1.00

0.00

1.00

2.00

3.00

4.00

5.00

6.00

7.00

8.00

-1.00 0.00 1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00

Experimental Densification (% TD)

Predicted Densification (% TD)

Banks

Freshley

Small

Freshley 1976

Figure 2.6-1 Comparison of densification predictions generated by the FUDENS model 

and experimental densification measurements determined in reactor. 

This data comparison demonstrates that basing densification on the sintering temperature 

provides a large degree of uncertainty. However, it should be noted that the Small (1987) data, 

which exhibited very high densification at very low burnup (less than 0.07 MWd/kgU), possessed 

a non-prototypic pore size distribution that was heavily weighted toward smaller diameters. 

Three of the four fuel types had 100 percent of their porosity in pores smaller than 2 m diameter, 

which is atypical of current fuel fabrication. 

Since in-reactor densification is strongly influenced by fuel fabrication, ex-reactor re-sintering 

tests are the preferred method for predicting in-reactor fuel densification (NRC, 1978). 

Ex-reactor resintering tests conducted by Freshly et al. (1976 and 1978) at 1600 and 1700°C 

provided reasonable estimate of in-reactor densification tests with standard deviations of 0.5 to 

0.6 percent of the theoretical density. 

2.6.2 Applicability and Uncertainty 

The densification correlation used in the subroutine FUDENS is applicable to fuels with pore size 

distributions that are similar to those included in the Freshley et al. (1976) study. Due to the 

scatter in the experimental data, it is difficult to establish a meaningful measure of uncertainty. 

2.6.3 References 

Rolstad, E., et al. 1974. “Measurements of the Length Changes of UO2 Fuel Pellets During 

Irradiation,” Enlarged HPG Meeting on Computer Control and Fuel Research, June 4-7, 1974. 
2.31 

NRC – U.S. Nuclear Regulatory Commission. 1978. An Acceptable Model and Related 

Statistical Methods for the Analysis of Fuel Densification. Regulatory Guide 1.126, Rev.1, U.S. 

Nuclear Regulatory Commission, Washington, D.C. 

Siefken, L.J., E.W. Coryell, E.A. Harvego, J.K. Hohorst. 2001. SCDAP/RELAP5/MOD 3.3 Code 

Manual MATPRO-A Library of Materials Properties for Light-Water-Reactor-Accident Analysis, 

NUREG/CR-6150, Vol. 4, Rev. 2, INEL-96/0422, Idaho National Engineering and 

Environmental Laboratory, Idaho Falls, ID. 

Hamming, R.W. 1971. Introduction to Applied Numerical Analysis, New York: McGraw-Hill 

Book Company, Inc. 

Freshley, M.D., D.W. Brite, J.L. Daniel, and P.E. Hart. 1976. “Irradiation-Induced Densification 

of UO2 Pellet Fuel,” Journal of Nuclear Materials, 62:138-166. 

Banks, D.A. 1974. “Some Observations of Density and Porosity Changes in UO2 Fuel Irradiated 

in Water-Cooled Reactors,” Journal of Nuclear Materials, 54:97-107. 

Freshley, M.D., D.W. Brite, J.L. Daniel, and P.E. Hart. 1978. “Irradiation-Induced Densification 

and PuO2 Particle Behavior in Mixed-Oxide Pellet Fuel,” Journal of Nuclear Materials, 81:63-

92. 

Small, G.J. 1987. “Densification of Uranium Dioxide at Low Burn-Up,” Journal of Nuclear 

Materials, 148:302-315. 

2.7 Fuel Swelling (FSWELL) 

The subroutine FSWELL calculates fuel swelling, which is caused by the buildup of solid and 

gaseous fission products during irradiation. Fuel swelling (FSWELL) is combined with 

creep-induced elongation (FCREEP), densification due to pressure sintering (FHOTPS), and 

irradiation (FUDENS) to calculate the overall dimensional changes in fuel. 

The gaseous swelling correlation in FSWELL is not used in FRAPCON-3 because it significantly 

overpredicts swelling and its effects on cladding deformation. The FSWELL correlation for solid 

swelling is used in both MATRPO and FRAPCON-3. However, the FSWELL solid swelling 

correlation used in FRAPCON-3.4 has evolved from the original model described in MATPRO. 

The FSWELL subroutine is not used in FRAPTRAN 1.4 because FRAPTRAN 1.4 is intended for 

shorter time scales than those that result in significant solid swelling. 

2.7.1 Model Development and Comparisons 

In MATPRO, the correlation used by FSWELL to calculate swelling due to the buildup of solid 

fission products is give by Equation 2.7.1. 

(2.7-1) 

where 

SS = fractional volume change due to solid fission products 

(m3

 volume change/m3 fuel) 

S

29

S B10x5.2S  
2.32 

BS = burnup during a time step (fissions/m3

) 

At temperatures below 2800K, swelling due to the buildup of gaseous fission products is given by 

Equation 2.7-2. 

 

(2.7-2) 

where 

Sg = fractional volume change due to gas fission products (m3

 volume change/m3 fuel) 

T = temperature (K) 

B = total burnup of fuel (fissions/m3

) 

For temperatures greater than 2800K, Sg is zero because the gas that causes swelling is assumed 

to have released. The correlations presented in Equations 2.7-1 and 2.7-2 are used to determine 

swelling due to solid and gaseous fission product buildup, respectively, in the FSWELL subcode 

described by MATPRO. 

Unlike the MATPRO version of FSWELL, FRAPCON-3.4 does not include a correlation for 

swelling due to the buildup of gaseous fission products. Results of an NRC-sponsored test rod 

were compared with predictions made by FRAPCON-3 and revealed that swelling models that 

included gaseous fission product buildup tend to overpredict test results at burnup levels up to 

10 GWd/MTU (Lanning and Bradley, 1984; Lanning et al., 1997). It was determined that 

swelling could be modeled better by eliminating the contribution from gaseous fission products 

and increasing the contribution of solid fission products from 0.699 percent (MATPRO and 

FRAPCON-2) to 0.77 percent (FRAPCON-3) per 10 GWd/MTU. 

In FRAPCON-3.4, the swelling due to solid fission products was decreased to 0.62 percent per 

10 GWd/MTU for burnup less than 80 GWd/MTU and increased to 0.86 percent per 

10 GWd/MTU for burnup greater than 80 GWd/MTU. The correlations used to calculate fuel 

swelling under these conditions are given in Equations 2.7-3 and 2.7-4, respectively. 

For burnup < 80 GWd/MTU 

(2.7-3) 

For burnup > 80 GWd/MTU 

(2.7-4) 

where 

soldsw = fractional volume change due to solid fission products (m3

 volume change/m3

fuel) 

sigswell = user-defined parameter that incrementally changes error 

bus = fuel burnup during time step (analogous to BS in MATPRO model) 

fdens = initial pellet density (kg/m3

) 

bu = burnup at end of time step (MWs/kgU) 

bul = burnup at end of previous time step (MWs/kgU) 

and 

       

S

56 73.11 B10x0.8T28000162.0

g eT280010x8.8S e B 27     

soldsw 10x315.2(*bus )10x315.2*sigswell 23 24  

soldsw 10x211.3(*bus )10x211.3*sigswell 23 24  
2.33 

The variable “bus” in Equations 2.7-3 and 2.7-4 is analogous to the variable “BS” in Equation 

2.7-1. However, “bus” permits the fractional volume change due to solid fission products to be 

determined as a function of burnup in MWs/kgU as opposed to fissions/m3

. 

A comparison between the MATPRO and the FRAPCON-3.4 versions of FSWELL is presented 

in Figure 2.7-1 along with experimental data from immersion density measurements. Based on 

this figure, it can be seen that the FRAPCON-3.4 model predicts less swelling than the MATPRO 

model, which tends to overpredict the experimentally determined data points. 

0.00

0.50

1.00

1.50

2.00

2.50

3.00

3.50

4.00

4.50

0 10 20 30 40 50 60 70 80

Burnup (GWd/MTU)

 V/V (%)

MATPRO FRAPCON-3.4 Garde (1986)

Newman (1986) Smith (1994) Dideon and Bain (1983)

Figure 2.7-1 Data comparison between MATPRO and FRAPCON-3.4 versions of 

FSWELL. 

For burnup > 80 GWd/MTU, a higher swelling rate is used in the FRAPCON-3.4 version of 

FSWELL. Data comparisons between FRAPCON-3.4 predictions and data from instrumented 

fuel assemblies tested at Halden are presented in Figure 2.7-2 for varying levels of burnup. These 

Halden swelling data are based on measurements of fuel stack length and cladding diameter 

change with burnup and converted to volume swelling assuming isotropic swelling. The 

diametral cladding change data versus burnup is particularly valuable because it implicitly 

includes the effects of cladding restraint on fuel swelling. 

 bulbu*10x974.2*fdensbus  10  
2.34 

0

0.2

0.4

0.6

0.8

1

1.2

1.4

IFA-633IFA-610

IFA-534.14

IFA-629.4IFA-629.3IFA-597.3

IFA-515.10IFA-515.10

IFA 519.9

IFA-504

IFA-655 Rods 1

&3

IFA-655 Rods 2

&4

Swelling Rate, %DV/V per 10 GWd/MTM

99 GWd/MTU99 GWd/MTU17 GWd/MTU- 62 GWd/MTU- 68 GWd/MTU- 65 GWd/MTU- 62 GWd/MTU- 90 GWd/MTU- 75 GWd/MTU- 78 GWd/MTU- 90 GWd/MTU- 70 GWd/MTUFigure 2.7-2 Comparison between swelling rates determined for various instrumented fuel 

assemblies tests at Halden and swelling rates predicted by FRAPCON-3.4. 

2.7.2 Applicability and Uncertainty 

Based on the data from which it is derived, the correlation used in FSWELL is applicable to fuels 

with burnups up to 100 GWd/MTU. Below 80 GWd/MTU, a 1-sigma uncertainty of 0.08 percent 

V/V per 10 GWd/MTU has been calculated. Above 80 GWd/MTU, a 1 sigma uncertainty of 

0.16 percent V/V per 10 GWd/MTU has been calculated. 

2.7.3 References 

Lanning, D.D., and E.R. Bradley. 1984. Irradiation History and Interim Postirradiation Data 

for IFA-432, NUREG/CR-3071, Pacific Northwest Laboratory for the U.S. Nuclear Regulatory 

Commission, Washington, D.C. 

Lanning, D.D. C.E. Beyer, C.L. Painter. 1997. FRAPCON-3: Modifications to Fuel Rod 

Material Properties and Performance Models for High-Burnup Application, NUREG/CR-6534, 

Vol. 1, Pacific Northwest National Laboratory for the U.S. Nuclear Regulatory Commission, 

Washington, D.C. 

Garde, A.M. 1986. Hot Cell Examination of Extended Burnup Fuel from Fort Calhoun, 

DOE/ET-34030-11 (CEND-427), Combustion Engineering, Windsor, CT. 

Newman, L.W. 1986. The Hot Cell Examination of Oconee 1 Fuel Rods after Five Cycles of 

Irradiation, DOE/ET-34212-36, Babcock and Wilcox, Lynchburg, VA. 
2.35 

Smith, G.P., R.C. Pirek, H.R. Freeburn, and D. Schrire. 1994. The Evaluation and 

Demonstration of Methods for Improved Nuclear Fuel Utilization, DOE/ET-34-013-15, 

Combustion Engineering, Windsor, CT. 

Dideon, C.G., and G.M. Bain. 1983. Fuel Performance under Extended Burnup for the B&W 

15x15 Design, DOE/ET-34212-38, Babcock and Wilcox, Lynchburg, VA. 
3.1 

3.0 Cladding Material Properties 

Material property correlations for fuel claddings are described in the following subsections. 

Unless otherwise specified, the correlations below are applicable to Zircaloy-2, -4, ZIRLO, and 

M5. In addition, several FRAPTRAN 1.4 subroutines include additional correlations to describe 

a Russian zirconium alloy containing 1 wt% niobium (ZrNb-1). The subroutine used to 

determine the material correlation in either the FRAPCON-3 or FRAPTRAN codes is given in 

parentheses. 

3.1 Cladding Specific Heat (CCP) 

The specific heat subcode, CCP, determines the true specific heat at constant pressures for 

cladding. Specific heat calculations are based on interpolation of measured data. The correlation 

in FRAPCON-3.4, FRAPTRAN 1.4, and MATPRO is applicable to Zircaloy-2, Zircaloy-4, 

ZIRLO, and M5 alloys. However, FRAPTRAN 1.4 contains correlations to calculate the specific 

heat of ZrNb-1. This calculation can be based on data collected at either fast or slow heating 

rates and the user can designate the appropriate data set. 

3.1.1 Model Development and Comparisons 

The CCP subcode requires temperature as an input to calculate specific heat. For the alpha phase 

of the Zircaloy alloys (temperature less than 1090K), CCP returns linear interpolations for the 

points listed in Table 3.1-1. These data points are based on precise data taken by Brooks and 

Stansbury (1966) with a Zircaloy-2 sample that had been vacuum-annealed at 1075K to remove 

hydrogen, which would have otherwise affected the measurement. 

The standard errors associated with this interpolation technique differ between MATPRO and 

FRAPCON-3.4/FRAPTRAN 1.4. In MATPRO, the standard error of the CCP interpolation was 

based on 90 points in the Brooks and Stansbury (1966) database and was found to be temperature 

dependent. For the 57 data points between 300 and 800K, the standard error is 1.1 J/kg*K. 

Between 800 and 1090K, it is 2.8 J/kg*K. The FRAPCON-3.4/FRAPTRAN 1.4 codes assume a 

standard error of 10 J/kg*K below 1090K. 

For temperatures from 1090 to 1300K (where Brooks and Stansbury do not report results), values 

of specific heat proposed by Deem and Eldridge (1967) are adopted. The Deem and Eldridge 

(1967) values are based on measurements of enthalpy and temperature which provide 

considerably less precise specific heat data than the results of Brooks and Stansbury (1966). 

As a result, the MATPRO standard error estimated from the Deem and Eldridge (1967) data in 

the region of 1090 through 1310K is 10.7 J/kg*K. This standard error is a measure only of the 

precision of the fit, since only a single data source is used. The standard error in 

FRAPCON-3.4/FRAPTRAN is assumed to be 25 J/kg*K between 1090 and 1300K. 

Above the alpha + beta to beta transformation temperature (about 1250K) and up to about 1320K, 

a constant value of 355.7 J/kg*K was reported by Deem and Eldridge (1967). This value agrees 

well with a value of 365.3 reported by Coughlin and King (1950) for pure beta zirconium. The 

standard error of specific heat calculations made above 1300K is assumed to be 100 J/kg*K in the 

FRAPCON-3.4/FRAPTRAN 1.4 codes. 
3.2 

In addition to Zircaloy alloys, FRAPTRAN 1.4 has been modified to include specific heat 

calculations for ZrNb-1. These calculations are based on data collected at two different heating 

rates. Depending on the user input, specific heat calculations can be based on either the fast or 

the slow heating rate data. The data used to interpolate specific heat values for the ZrNb-1 alloy 

are presented in Table 3.1-2. There is no standard error described for specific heat calculations 

based on these data sets. 

Table 3.1-1 Specific heat capacity database for Zircaloy-2, Zircaloy-4, ZIRLO, and M5. 

 Standard Error 

Temperature

(K)

Specific Heat 

Capacity

(J/kg*K)

Source MATPRO

(J/kg*K)

FRAPCON-3 / 

FRAPTRAN

(J/kg*K)

----------------------------Alpha Phase---------------------------- 

300 281 Brooks and Stansbury 1.1 10 

400 302 Brooks and Stansbury 1.1 10 

640 331 Brooks and Stansbury 1.1 10 

1090 375 Brooks and Stansbury 2.8 10 

----------------------------Beta Phase---------------------------- 

1093 502 Deem and Eldridge 10.7 25 

1113 590 Deem and Eldridge 10.7 25 

1133 615 Deem and Eldridge 10.7 25 

1153 719 Deem and Eldridge 10.7 25 

1173 816 Deem and Eldridge 10.7 25 

1193 770 Deem and Eldridge 10.7 25 

1213 619 Deem and Eldridge 10.7 25 

1233 469 Deem and Eldridge 10.7 25 

1248 356 Deem and Eldridge 10.7 25 

2098 356 Coughlin and King 100 100 

2099 356 Coughlin and King 100 100 

Table 3.1-2 Specific heat capacity database for ZrNb-1 alloys. 

 

Temperature

(K)

Specific Heat Capacity

(J/kg*K)

Source

--------- Low Heating Rate ( < 1000 K/s) --------- 

280 345 Volkov B. Yu et al. 

473 360 Volkov B. Yu et al. 

573 370 Volkov B. Yu et al. 

673 380 Volkov B. Yu et al. 

773 383 Volkov B. Yu et al. 

873 385 Volkov B. Yu et al. 

883 448 Volkov B. Yu et al. 

973 680 Volkov B. Yu et al. 

1025 816 Volkov B. Yu et al. 

1073 770 Volkov B. Yu et al. 
3.3 

1153 400 Volkov B. Yu et al. 

1173 392 Volkov B. Yu et al. 

1248 356 Volkov B. Yu et al. 

--------- High Heating Rate ( > 1000 K/s) --------- 

1100 412 Ljusternik V.E. et al. 

1110 420 Ljusternik V.E. et al. 

1120 480 Ljusternik V.E. et al. 

1134 600 Ljusternik V.E. et al. 

1142 1000 Ljusternik V.E. et al. 

1150 1400 Ljusternik V.E. et al. 

1155 1600 Ljusternik V.E. et al. 

1161 1400 Ljusternik V.E. et al. 

3.1.2 Applicability and Uncertainty 

The CCP subroutine is valid over the temperature range from which the data sets have been 

collected on zirconium alloys (300 to 2099K). This correlation is applied to Zircaloy-2, -4, 

ZIRLO, and M5. However the standard error increases with increasing temperature as indicated 

in Table 3.1-1. 

3.1.3 References 

Brooks, C.R., and E.E. Stansbury. 1966. “The Specific Heat of Zircaloy-2 from 50 to 700°C,” 

Journal of Nuclear Materials, 18:223. 

Deem, H.W., and E.A. Eldridge. 1967. Specific Heats of Transformation of Zirclaoy-2 and Low 

Nickel Zircaloy-2, USAEC BM1-1803, Battelle Memorial Institute, Columbus, OH. 

Coughlin, J.P., and E.G. King. 1950. “High-Temperature Heat Contents of Some ZirconiumContaining Substances,” Journal for the American Chemical Society, 72:2262. 

Volkov, B. Yu., et al. 1989. Library of Subprograms on Physical and Mechanical Properties of 

the N1-Alloy Fuel Rod Cladding Material, Preprint IAE-4941/11, Institute of Automation and 

Electrometry, Moscow, Russia. 

Ljusternik, V.E., et al. 1993. “Experimental research of zirconium reactor materials thermal 

properties: 110-alloy,” M.J. Higher Temperature Thermal Physics, 31(4). 

3.2 Cladding Thermal Conductivity (CTHCON) 

The subroutine CTHCON is used to calculate cladding thermal conductivity, which is required 

for accurate predictions of fuel temperature. The thermal conductivity of the cladding is 

primarily a function of temperature. Other characteristics, such as residual stress levels, crystal 

orientation, and minor composition differences, may have secondary effects on thermal 

conductivity. The correlation used in CTCHCON to calculate cladding thermal conductivity is 

the same in MATPRO, FRAPCON-3.4, and FRAPTRAN 1.4. This correlation is applied to 

Zircaloy-2, -4, ZIRLO, and M5. In addition, FRAPTRAN 1.4 includes a correlation to calculate 

the cladding thermal conductivity of ZrNb-1. 
3.4 

3.2.1 Model Development and Comparisons 

Considering only temperature as the defining parameter, the thermal conductivity of Zircaloy for 

temperatures less than 2098K is described by CTHCON in Equation 3.2-1. The uncertainty of 

the conductivity calculations is given in Equation 3.2-2. 

(3.2-1) 

(3.2-2) 

For temperatures greater than or equal to 2098K, the thermal conductivity and uncertainty are 

given in Equations 3.2-3 and 3.2-4, respectively. 

(3.2-3) 

(3.2-4) 

where 

k = thermal conductivity of Zircaloy (W/m* K) 

T = temperature (K) 

k = standard deviation (W/m*K) 

Equation 3.2-1 predicts k very well from room temperature to the data limit of about 1800K 

(Figure 3.2-1) and may be extrapolated with some confidence to the melting point. The standard 

deviation (k) of the data with respect to this correlation appears to be temperature independent 

over the data range. 

The standard deviation of the constants in Equation 3.2-1 is 20 to 30 percent of the value of the 

constant. Jensen (1969) performed a parametric analysis of several variables involved in the 

estimation of fuel and cladding temperatures. Both steady-state and transient analysis showed 

that variations of  20 percent resulted in calculated cladding temperature variations of about 

2.8K. Fuel centerline temperatures are more sensitive to cladding thermal conductivity and 

showed variations of 28K. Similar findings were reported by Korber and Unger (1974). 

For temperatures less than 2133K, the thermal conductivity of the ZrNb-1 alloy is calculated by 

FRAPTRAN 1.4 with the correlation presented in Equation 3.2-5. 

(3.2-5) 

For temperatures above 2133K, the thermal conductivity of the ZrNb-1 cladding is assumed to be 

36 W/m*K. 

Figure 3.2-1 compares the CTHCON model as incorporated in MATPRO, FRAPCON-3.4, and 

FRAPTRAN-1.4 with the data on which it is based, extrapolated to the melting point of Zircaloy 

(2098K). 

Figure 3.2-2 compares the CTHCON model with the Zircaloy thermal conductivity and 

2 25 39 T10x67.7T10x45.1T10x09.251.7k      

01.1 k 

 36k

5 k 

000461843.0 T*  e0636.15k
3.5 

diffusivity data that were not included in development of the model. Thermal diffusivity data 

were converted to thermal conductivity using the recommended models in Fink (2000) for 

specific heat and density, respectively. Although different alloys were examined, there were no 

discernable differences between the alloys. The scatter in the data due to differences in alloy, 

form, or measurement direction is comparable to the variation between individual data sets. The 

scatter in the data is distributed relatively uniformly around the curve representing the CTHCON 

model, except for anomalously high Peggs data below 1000K and the anomalously low Bunnell 

et al. (1983) data above 1300K. 

0

10

20

30

40

50

60

70

0 500 1000 1500 2000 2500

Temperature (K)

Thermal Conductivity (W/m-K) 

MATPRO

Anderson et al. (1962)

Chirigos et al. (1961)

Feith (1966)

Lucks and Deem (1958)

Powers (1961)

Scott (1965)

Figure 3.2-1 CTHCON model and data used for development. 
3.6 

0

10

20

30

40

50

60

70

0 500 1000 1500 2000 2500

Temperature (K)

Thermal Conductivity (W/m-K) 

MATPRO

Krett and Cleveland (1997) - China

Krett and Cleveland (1997) - AECL

Gilchrist (1976)

Bunnell et al. (1983)

Murabayashi et al. (1975)

Peggs et al. (1976)

Maglic et al. (1994)

Figure 3.2-2 CTHCON model and data collected since development. 

3.2.2 Applicability and Uncertainty 

The correlation used in CTHCON to calculate the cladding thermal conductivity may be applied 

to temperatures up to the melting temperature of Zircaloy (2098K). The uncertainty through this 

range of temperature is given by Equation 3.2-2. Above the melting temperature, the uncertainty 

increases to the value given by Equation 3.2-4. 

3.2.3 References 

Jensen, S.E. 1969. Parametric Studies of Fuel Pin Temperature Response, IDO-17295, Philips 

Petroleum Company, Atomic Energy Division, Idaho Falls, ID. 

Korber, H., and H.E. Unger. 1974. “Sensitivity Study on Core Heatup and Meltdown by 

Variation of Heat Conductivity and Thermal Emissivity,” Transactions for the American Nuclear 

Society, 18:234-235. 

Fink, J.K. 2000. “Zircaloy Thermal Conductivity,” International Nuclear Safety Center Material 

Properties Database, http://www.insc.anl.gov/matprop/zircaloy/index.php, Argonne National 

Laboratory, Argonne, IL. 

Bunnell, L.R., J.L. Bates, and G.B. Mellinger. 1983. “Some High-Temperature Properties of 

Zircaloy-Oxygen Alloys,” Journal of Nuclear Materials, 116(2-3):219-232. 

Anderson, W.K., C.J. Beck, A.R. Kephar, and J.S. Theilacker. 1962. “Zirconium Alloys,” 

Reactor Structural Materials: Engineering Properties as Affected by Nuclear Reactor Service, 
3.7 

ASTM-STP-314, pp. 62-93, American Society for Testing and Materials, West Conshohocken, 

PA. 

Chirigos, J.N., et al. 1961. “Development of Zircaloy-4,” Fuel Element Fabrication, New York: 

Academic Press, pp. 19-55. 

Feith, A.D. 1966. Thermal Conductivity and Electrical Resistivity of Zircaloy-4, GEMP-669, 

October 1966, General Electric Company, Missile and Space Division, Cincinnati, OH. 

Lucks, C.F., and H.W. Deem. 1958. Progress Relating to Civilian Applications during June, 

1958, R.W. Dayton and C.R. Tipton, Jr., eds., BMI-1273, pp. 7-9, Battelle Memorial Institute, 

Columbus, OH. 

Powers, A.E. 1961. Application of the Ewing Equation for Calculating Thermal Conductivity 

from Electrical Conductivity from Electrical Conductivity, KAPL-2146, Knolls Atomic Power 

Laboratory, General Electric Company, Schenectady, NY. 

Scott, D.B. 1965. Physical and Mechanical Properties of Zircaloy 2 and 4, WCAP-3269-41, 

pp. 5, 9, Westinghouse Electric Corporation, Atomic Power Division, Pittsburgh, PA. 

Krett, V., and J. Cleveland. 1997. Thermophysical Properties of Materials for Water Cooled 

Reactors, IAEA-TECDOC-949, pp. 67-84, International Atomic Energy Agency, Vienna, 

Austria. 

Gilchrist, K.E. 1976. “Thermal Property Measurements on Zircaloy-2 and Associated Oxide 

Layers up to 1200°C,” Journal of Nuclear Materials, 62:257-264. 

Murabayashi, M., S. Tanaka, and Y. Takahashi. 1975. “Thermal Conductivity and Heat 

Capacity of Zircaloy-2, -4 and Unalloyed Zirconium,” Journal of Nuclear Science and 

Technology, 12(10):661-662. 

Peggs, I.D., A.M. Stadnyk, and D.P. Godin. 1976. “Thermophysical Properties of ZirconiumAlloy Fuel Channel Components,” High Temperatures-High Pressures, 8(4):441-450. 

Maglic, K.D., N.L.j. Perovic, and A.M. Stanimirovic. 1994. “Calorimetric and Transport 

Properties of Zircalloy 2, Zircalloy 4, and Inconel 625,” International Journal of Thermophysics, 

15(4):741-755. 

3.3 Cladding Oxide Thermal Conductivity (ZOTCON) 

The ZOTCON subroutine calculates the thermal conductivity of the zirconium oxide layer that 

forms on zirconium alloys. Cladding temperature is the only parameter used to calculate 

zirconium oxide thermal conductivity. The correlation used in the ZOTCON subroutine differs 

between MATPRO and the FRAPCON-3.4/FRAPTRAN 1.4 codes. The MATPRO correlation is 

based on an arbitrary correction to an earlier correlation that was fit to a large database of oxide 

thermal conductivity measurements. This database contained measurements from a variety of 

zirconium oxide materials and test methods. The correction was applied to better fit data that was 

more representative of the oxide layer that forms on zirconium alloys. The correlation used in 

FRAPCON-3.4 and FRAPTRAN 1.4 is based on measurements obtained by Kingery (1954) from 

fully dense and porous (87 percent TD) zirconium oxides. 
3.8 

3.3.1 Model Development and Comparisons 

Oxide thermal conductivity models in MATPRO and FRAPCON-3.4/FRAPTRAN 1.4 are only a 

function of temperature. After the cladding temperature is established, the thermal conductivity 

of the oxide layer can be calculated using the correlation determined from either the MATPRO or 

the FRAPCON-3.4/FRAPTRAN 1.4 databases. 

The MATPRO (2001) correlation is derived from a database that includes data collected by 

Adams (1954), Maki (1973), Lapshov and Bashkatov (1973), and Gilchrist (1976). The 

correlation used in FRAPCON-3.4/FRAPTRAN 1.4 is derived from data collected by Kingery 

(1954). The FRAPCON-3.4/FRAPTRAN correlation and database was adopted from an earlier 

version of MATPRO (Hagrman et al., 1981). 

The database for the MATPRO (2001) model includes data from a variety of zirconium oxide 

materials and test conditions. The Maki (1973) data include two samples oxidized in steam over 

a small temperature range and show a sharp increase in conductivity between 400 and 500K. The 

principal recommendation for the data is that they were taken with black oxide from Zircaloy 

tubes. Two sets of data attributed to Waldman by Maki are also shown in the table. The data of 

Lapshov and Bashkatov (1973) are from films formed by plasma sputtering of zirconium dioxide 

on tungsten substrates. However, these data may not be representative of Zircaloy cladding oxide 

thermal conductivity because sputtered coatings are quite porous. In addition, these films do not 

have the same O/M ratio as cladding oxide and may not adhere well to the substrate. The 

Gilchrist (1976) data includes two types of oxide films. One is a nodular oxide and the other a 

black oxide characteristic of the kinds of layers usually reported in high-temperature tests with 

cladding. The nodular oxide thermal conductivities are much lower than the black oxide thermal 

conductivities, and both kinds of oxide have conductivities that are significantly lower than the 

stabilized Zircaloy dioxide conductivities reported by Adams (1954). Considerable uncertainty is 

reported by Gilchrist (1976) because of difficulty in measuring oxide film thickness. 

The database for the FRAPCON-3.4 and FRAPTRAN 1.4 model is based on data from Kingery 

(1954). The data are from two different samples; one was a bulk sample of 100 percent dense 

ZrO2 and the other with approximately 13 percent porosity. 

Each database was fit with a different model to relate the thermal conductivity of the cladding 

oxide layer to temperature. The equations used to describe these models are presented in 

Equations 3.3-1 and 3.3-2 for MATPRO and FRAPCON-3.4/FRAPTRAN 1.4 (i.e., FRAP), 

respectively. These models relate the thermal conductivity of the cladding oxide layer, K 

(W/m*K), to cladding temperature, T (K). 

(3.3-1) 

(3.3-2) 

Both databases are presented in Figure 3.3-1 along with both cladding oxide thermal conductivity 

models. In addition, a more recent data set from Gilchrist (1979) is included. Figure 3.3-1 

reveals that the FRAP model fits the Kingery data from which it was derived and agrees well with 

the measurements performed by Adams (1954). The FRAP model tends to overestimate both the 

Gilchrist data sets and the data from Lapshov and Bashkatov (1973). However, as stated earlier, 

the earlier Gilchrist (1976) data is associated with considerable uncertainty and the data from 

K T10x81.1835.0 4

MATPRO

 

))10x946.1*T10x43.6(*T10x41.2(*T9599.1K 4 7 10

FRAP

     
3.9 

Lapshov and Bashkatov (1973) is likely to be low due to different O/M ratios and high porosity in 

the sputtered films. 

The MATPRO model was initially fit to the data collected by Adams (1954), but it was adjusted 

to the form presented in Equation 3.3-1 to better fit the Gilchrist (1976) data for black oxide. The 

rationale for this correction was that the black oxide was more prototypic of the cladding oxide 

layer and was achieved by dividing the original equation by two. 
3.10 

Figure 3.3-1 Oxide thermal conductivity measurements as a function of temperature compared with the MATPRO a

3.4/FRAPTRAN 1.4 (FRAP) cladding oxide thermal conductivity models. 

0123456789

300 500 700 900 1100 1300 1500 1700 

Temperature (K)

Thermal Conductivity (W/m*K)

Adams Maki Maki (Waldman data)Lapshov and Bashkatov Gilchrist (black oxide) Gilchrist (nodular oxide)Kingery (13% Porosity) Kingery (100% Dense) Gilchrist (1979)

MA
3.11 

3.3.2 Applicability and Uncertainty 

The subroutine ZOTCON is applicable over the temperature range for which the data was 

collected, 375 to 1673K. Since the relatively high density of the bulk ZrO2 materials bias the 

thermal conductivity measurements high, it is difficult to establish a meaningful measurement of 

uncertainty. 

The MATPRO model is applicable from 300K to the melting temperature (2973K). Upon 

melting the thermal conductivity, KLiquid (W/m*K), is assumed to be described by Equation 3.3-3 

(Hagrman et al., 1981). 

(3.3-3) 

The standard error for the MATPRO model is approximately 0.75 W/m*K, but is assumed to be 

within 10 percent for materials that are known to be dense (>87%TD) ZrO2. 

3.3.3 References 

Adams, M. 1954. “Thermal Conductivity: III, Prolate Spheroidal Envelope Method,” Journal of 

the American Ceramic Society, 37:74-79. 

Maki, H. 1973. “Heat Transfer Characteristics of Zircaloy-2 Oxide Film,” Journal of Nuclear 

Science and Technology, 10:107-175. 

Lapshov, V. N., and A. V. Bashkatov. 1973. “Thermal Conductivity of Coatings of Zirconium 

Dioxide Applied by the Plasma Sputtering Method,” Heat Transfer-Soviet Research, 5:19-22. 

Gilchrist, K. E. 1976. “Thermal Property Measurements on Zircaloy-2 and Associated Oxide 

Layers,” Journal of Nuclear Materials, 62:257-264. 

Kingery, W.D., et al. 1954. “Thermal Conductivity: X, Data for Several Pure Oxide Materials 

Corrected to Zero Porosity,” Journal of the American Ceramic Society, 37:107-100. 

Hagrman, D.L., et al. 1981. MATPRO-Version 11 (Revision 2). A Handbook of Materials 

Properties for Use in the Analysis of Light Water Reactor Fuel Rod Behavior, NUREG/CR-0479 

(TREE-1280, Rev. 2), EG&G Idaho, Inc., Idaho Falls, ID. 

Gilchrist, K.E. 1979. “Thermal Conductivity of Oxide Deposited on Zircaloy Fuel Tube 

Material – A Continuation of Previous Work,” Journal of Nuclear Materials, 82:193-194. 

3.4 Cladding Surface Emissivity (ZOEMIS) 

The subcode ZOEMIS returns the cladding surface emissivity, which is directly proportional to 

the radiant heat transfer from the cladding surface during an abnormal transient. The ZOEMIS 

model described in MATPRO is the same as the model used in FRAPCON-3.4 and FRAPTRAN 

1.4. 

Liquid  4.1K
3.12 

3.4.1 Model Development and Comparisons 

Surface emissivities are significantly affected by surface layers on the cladding. For cladding 

with thin oxide coatings, the oxide surface thickness is only a few wavelengths of near infrared 

radiation and is partly transparent. Oxide thickness is an important parameter for these thin 

coatings. Thicker oxide layers are opaque, so the oxide thickness is not as important as the nature 

of the outer oxide surface, which is affected by temperature and by chemical environment. The 

effect of temperature has been modeled, but variations in crud on the external cladding surface 

and chemical reaction products on the inside surface are not modeled explicitly. 

The model for emissivity was constructed by considering measured emissivities reported by 

several investigators. Expressions used to predict the emissivity of Zircaloy cladding surfaces are 

summarized below. 

When the cladding surface temperature has not exceeded 1500K, emissivities are modeled by 

Equations 3.4-1 and 3.4-2. Equation 3.4-1 is used for oxide layers less than 3.88x10-6 m thick 

and Equation 3.4-2 is used for oxide layers equal to or greater than 3.88x10-6 m thick. Both 

equations relate the hemispherical emissivity, 1 (unitless), to the oxide layer thickness, d (m). 

 (3.4-1) 

 (3.4-2) 

When the maximum cladding temperature has exceeded 1500K, emissivity is taken to be the 

larger of 0.325 and the result of Equation 3.4-3. Equation 3.4-3 relates the emissivity above 

1500K, 2 (unitless), to 1 and the maximum cladding temperature, T (K). 

 (3.4-3) 

Comparisons between data and the high and low temperature models are presented in Figure 

3.4-1. 

d10x1246.0325.0 6

1 

d0.50808642.0 1  





 



  

300

T1500

exp 12
3.13 

0.00.10.20.30.40.50.60.70.80.91.0

0 20 40 60 80 100 120 

Surface Layer Thickness (mm)

Emissivity (unitless) Murphy and Havelock Juenke and Sjodahl

Burgoyne and Garlick Uniform Oxide Burgoyne and Garlick Nodular OxidBurgoyne and Garlick with Crud Model Prediction < 1500 KModel Prediction > 1500K

Figure 3.4-1 Model-to-data comparison for cladding oxide emissivity. 
3.14 

3.4.2 Applicability and Uncertainty 

The absolute standard error, 1, expected from the use for Equation 3.4-1 to predict emissivity in 

a reactor when cladding surface temperature has never exceeded 1500K is given by Equation 

3.4-4. 

 (3.4-4) 

When the maximum cladding temperature exceeds 1500K, the expected standard error is 

estimated by 2, which is described by Equation 3.4-5. 

 (3.4-5) 

If Equations 3.4-3 and 3.4-5 predict values of 2 and 2 that fall inside the range of physically 

possible values of emissivity (0.0 – 1.0), the value 2 is returned as the expected standard error. 

If the prediction 2 + 2 is greater than 1 or if 2 - 2 is less than 0, the standard error of Equation 

3.4-5 is modified to limit 2 + 2 at 1 and/or 2 - 2 at 0. 

3.4.3 References 

Juenke, E.F., and L.H. Sjodahl. 1968. “Physical and Mechanical Properties: Emittance 

Measurements,” Advanced Pressure Vessel Materials, GEMP-1008, pp. 239-242, General 

Electric Co., Missile and Space Division, Cincinnati, OH. 

Murphy, E.V., and F. Havelock. 1976. “Emissivity of Zirconium Alloys in Air in the 

Temperature Range 100-400°C,” Journal of Nuclear Materials, 60:167-176. 

Burgoyne, T.B. and A. Garlick. 1976. “The effect of oxidation and crud deposition on the 

emissivity of Zircaloy-2 cladding.” Paper presented to the OECD-CSNI Meeting on the Behavior 

of Water Reactor Fuel Elements under Accident Conditions, Spåtnid, Norway, Sept. 13-16, 1976. 

3.5 Cladding Thermal Expansion (CTHEXP) 

The subroutine CTHEXP returns the axial and diametral components of thermal expansion in the 

cladding as a function of temperature. The MATPRO model calculates these strain components 

for single crystal Zircaloy and applies these to polycrystalline cladding materials through the use 

of pole figures collected over a cladding section. The FRAPCON-3.4/FRAPTRAN 1.4 models 

use the correlation from Mehan and Wiesinger (1961) for temperatures below 1073K and a 

constant value of thermal expansion above 1273K. Values of thermal expansion between 1073 

and 1273K are determined in the FRAPCON-3.4/FRAPTRAN 1.4 codes by linear interpolation. 

In addition, the FRAPTAN 1.4 model includes a calculation of the thermal expansion for ZrNb-1. 

3.5.1 Model Development and Comparisons 

In the MATPRO model, a total of six correlations that are functions of temperature only are used 

to find single crystal thermal strains. In addition, basal plane symmetry (11 = 22) is assumed. 

1.0 1  





 



  

300

1500T exp1.0 2
3.15 

The model was developed for as fabricated Zircaloy-4, but comparisons with Zircaloy-2 and 

zirconium data also show good agreement for these materials. 

The correlations for single crystal thermal strains are: 

For 300 < T < 1083K 

(3.5-1) 

(3.5-2) 

where 

11 = circumferential thermal expansion (m/m) 

33 = axial thermal expansion (m/m) 

For 1083  T < 1244K 

(3.5-3) 

(3.5-4) 

where the arguments of the cosines are in radians. 

For 1244  T < 2098K 

(3.5-5) 

(3.5-6) 

For temperatures  2098K, consideration of the volume change associated with melting is 

required. Since no data were found for the MATPRO model, a typical 2 percent volume increase 

at melt was assumed (Siefken et al., 2001). The expressions used for the thermal strain in liquid 

zirconium (temperatures  2098K) are thus: 

(3.5-7) 

p = thermal expansion strain in liquid Zircaloy (m/m) 

11 = circumferential thermal expansion strain of a single crystal of Zircaloy at 2098K 

(m/m) 

33 = axial thermal expansion strain of a single crystal of Zircaloy at 2098K (m/m) 

where the 11 and 33 are calculated by Equations 3.5-3 and 3.5-4, respectively. 

6 3

11 10x485.1T10x95.4    

5 3

33 10x78.3T10x26.1    

3

11 10x

161

1083T cos09822.177763.2 





 



 



 





   

3

33 10x

161

1083T cos09822.177763.2 





 



 



 





   

6 2

11 10x04.1T10x7.9    

6 3

33 10x4.4T10x76.9    

0067.0

3

1

3

2

p 11 33 
3.16 

To obtain cladding strains from these single crystal strains, it is necessary to do a volume 

weighted averaging of crystalline orientation over the entire cladding section. Such an averaging 

requires the collection of pole figures and is described in greater detail in the MATPRO (2001) 

handbook. The input parameters for typical cladding are not readily available. 

The model for cladding thermal expansion in the FRAPCON-3.4/FRAPTRAN (FRAP) codes is 

different from the MATPRO model but provides similar predictions of expansion (< 6 percent). 

The data for the FRAP correlation used from room temperature to 1273K was taken from Mehan 

and Wiesinger (1961), Scott (1965), and Kearns (1965). Above 1273K, the coefficient of thermal 

expansion is the constant value of 9.7x10-6, as recommended by Lustman and Kerze (1955). 

Between 1073 and 1273K (approximately the alpha-beta transition range for Zircaloy), the 

thermal expansion components are determined by linear interpolation. 

The correlations used to calculate the axial and diametral components of thermal expansion in the 

cladding between room temperature and 1073K are presented in Equations 3.5-8 and 3.5-9, 

respectively. Strain is given a function of temperature, T (°C). 

(3.5-8) 

(3.5-9) 

The correlations used to calculate the axial and diametral components of thermal expansion in the 

cladding above 1273K are presented in Equations 3.5-10 and 3.5-11, respectively. Strain is given 

a function of temperature, T (°C). 

(3.5-10) 

(3.5-11) 

Model-to-data comparisons for the MATPRO and FRAPCON-3.4/FRAPTRAN 1.4 models are 

presented in Figures 3.5-1 and 3.5-2 for axial and diametral thermal expansion, respectively. 

These comparisons illustrate that the axial and diametral components of cladding thermal 

expansions are being underpredicted. These correlations will be examined further and possibly 

revised in a future revision of FRAPCON-3.4/FRAPTRAN 1.4. 

axial T 5 6 104410.4105060.2    

diametral T 5 6 107210.6103730.2    

axial T 3 6 1070.910300.8    

diametral T 3 6 1070.910800.6    
3.17 

0.0E+00

1.0E-03

2.0E-03

3.0E-03

4.0E-03

5.0E-03

6.0E-03

7.0E-03

8.0E-03

9.0E-03

1.0E-02

200 400 600 800 1000 1200 1400 1600 1800

Temperature (K)

Thermal Expansion (m/m)

FRAP Axial

MATPRO Axial

Bunnel Axial

Kearns, Scott, and Mehan Axial

Figure 3.5-1 Model-to-data comparison for FRAPCON-3.4/FRAPTRAN 1.4 (FRAP) and 

MATPRO axial thermal expansion. Data from cladding tubes. 

0.0E+00

1.0E-03

2.0E-03

3.0E-03

4.0E-03

5.0E-03

6.0E-03

7.0E-03

8.0E-03

9.0E-03

1.0E-02

200 400 600 800 1000 1200 1400 1600 1800

Temperature (K)

Thermal Expansion (m/m)

FRAP Diametral

MATPRO Diametral

Bunnel Diametral

Kearns, Scott, and Mehan Diametral

Figure 3.5-2 Model-to-data comparison for FRAPCON-3.4/FRAPTRAN 1.4 (FRAP) and 

MATPRO diametral thermal expansion. Data from cladding tubes. 

3.5.2 Applicability and Uncertainty 

The subroutine CTHEXP is applicable over the temperature range for which the experimental 

data was collected, which included temperatures up to 1100K. The uncertainty of the axial and 

diametral thermal expansion correlations has not been quantified. 
3.18 

3.5.3 References 

Mehan, R.L., and F.W. Wiesinger. 1961. Mechanical Properties of Zircaloy-2, KAPL-2110, 

General Electric Company, Knolls Atomic Power Laboratory, Schenectady, NY. 

Siefken, L.J., E.W. Coryell, E.A. Harvego, and J.K. Hohorst. 2001. SCDAP/RELAP5/MOD 3.3 

Code Manual MATPRO-A Library of Materials Properties for Light-Water-Reactor-Accident 

Analysis, NUREG/CR-6150, Vol. 4, Rev. 2, INEL-96/0422, Idaho National Engineering and 

Environmental Laboratory, Idaho Falls, ID. 

Scott, P.B. 1965. Physical and Mechanical Properties of Zircaloy-2 and -4, WCAP-3269-41, 

Westinghouse Electric Corporation, Atomic Power Division, Pittsburgh, PA. 

Kearns, J.J. 1965. Thermal Expansion and Preferred Orientation in Zircaloy, WAPD-TM-472, 

pp. 17-18, Bettis Atomic Power Laboratory, Pittsburgh, PA. 

Lustman and F. Kerze. 1955. The Metallurgy of Zirconium, McGraw-Hill Book Company, New 

York, p. 355. 

Bunnell, L.R., et al. 1977. High Temperature Properties of Zircaloy Oxygen Alloys, EPRI NP524, Battelle Pacific Northwest Laboratory, Richland, WA. 

3.6 Cladding Elastic Modulus (CELMOD) and Shear Modulus 

(CSHEAR) 

Elastic moduli are required to relate stresses to strains. The elastic moduli are defined by the 

generalized form of Hooke’s law as elements of the fourth rank tensor that relates the second rank 

stress and strain tensors below the yield point. In practice, cladding is frequently assumed to be 

an isotropic material. In such a case, only two independent elastic moduli are needed to describe 

the relation between elastic stress and strain: the Young’s modulus and the shear modulus. 

The subcodes CELMOD and CSHEAR are used in the MATPRO and the 

FRAPCON-3.4/FRAPTRAN 1.4 (FRAP) codes to determine the Young’s modulus and the shear 

modulus. The tensor from which these moduli are derived is calculated by CELAST, which is 

included in MATPRO and determines the compliance matrix for isotropic cladding. However, 

CELAST is not used in the FRAP codes. The subcode CELAST would only be required if the 

Young’s modulus and shear modulus for an anisotropic cladding is desired and, as stated 

previously, the cladding is assumed to be isotropic. 

 

Elastic moduli are affected primarily by temperature and oxygen content. Fast neutron fluence, 

cold work, and texture effects are also included in the models described herein, but they are not as 

important as temperature and oxygen content for typical LWR fuel rod cladding. The models are 

based primarily on data published by Bunnell et al. (1977), Fisher and Renken (1964), Armstrong 

and Brown (1964), and Padel and Groff (1976), since these data include the best description of 

texture for the temperature range in which they were used. Data from several other sources are 

used to evaluate the expected standard error of the CELMOD and CSHEAR codes and to 

estimate the effect of fast neutron fluence (Shober, 1957; Whitmarsh, 1962). To calculate 

Zircaloy elastic moduli at temperatures greater than the melting temperature of Zircaloy (2098K), 

the moduli are set to zero (actually, 1.0x10-10 is used to avoid dividing by zero). 
3.19 

3.6.1 Model Development and Comparisons 

CELMOD

The CELMOD subcode used in MATPRO is the same equations written in a different formas the 

one used in FRAPCON-3.4/FRAPTRAN 1.4 (FRAP). The differences in form between the codes 

will be discussed in this subsection. 

MATPRO uses Equations 3.6-1 and 3.6-2 to describe the Young’s modulus in the alpha phase 

and the beta phase, respectively. 

(3.6-1) 

(3.6-2) 

where 

Y = Young’s modulus for Zircaloy-2 and -4 with random texture (Pa) 

T = cladding temperature (K) 

K1 = modification to account for the effect of oxidation (Pa) 

K2 = modification to account for the effect of cold work (Pa) 

K3 = modification to account for the effect of fast neutron fluence (unitless) 

In the alpha plus beta phase, MATPRO determines Y by linearly interpolating between values 

calculated at the alpha to alpha plus beta and the alpha plus beta to beta phase boundaries. 

The MATPRO expressions used to model the effects of oxidation, cold work, and fast neutron 

fluence are presented in Equations 3.6-3 through 3.6-5. 

(3.6-3) 

(3.6-4) 

(3.6-5) 

where 

 = average oxygen concentration minus oxygen concentration of as-received 

cladding (kg oxygen/kg Zircaloy). As-received oxygen concentrations are so 

small (0.0012 kg oxygen/kg Zircaloy) that the exact magnitude of the as-received 

concentration will not affect the correlation predictions. 

C = cold work (unitless ratio of areas). 

 = fast neutron fluence (n/m2

). 

The standard error of CELMOD is 6.4x109

 Pa. Although FRAPCON-3.4 is coded to include this 

error, this value is not used. This standard error is not included in FRAPTRAN 1.4. 

321

11 7    K/)KKT10x475.510x088.1(Y

T10x05.410x21.9Y 10 7  

  )T10x912.510x61.6(K  11 8

1

C10x6.2K 10

2 





 



 

  3 25 10 exp12.088.0K
3.20 

The CELMOD subcode included in the FRAP codes differs only in formy from the MATPRO 

version. For instance, the expression for Young’s modulus in the alpha phase is presented in 

Equation 3.6-6. 

(3.6-6) 

where 

celmod = Young’s modulus (Pa) 

ctemp = cladding temperature (K) 

deloxy = input average oxygen concentration excluding oxide layer (kg oxygen/kg 

Zircaloy)(hardwired to zero in FRAPCON-3.4/FRAPTRAN 1.4) 

cwkf = input effective cold work (unitless ratio of areas) 

In Equation 3.6-6, c1, c3, and c2 are expressions that account for oxygen content, cold work, and 

fast neutron fluence respectively. 

Equations 3.6-7 through 3.6-9 describe the constants c1, c2, and c3, respectively. 

(3.6-7) 

(3.6-8) 

(3.6-9) 

For neutron fluences greater than 1x1022, c2 is given by Equation 3.6-10. 

(3.6-10) 

where 

fnck = input effective fast fluence (n/m2

) 

Both FRAPCON-3.4 and FRAPTRAN 1.4 use the same correlations to determine Young’s 

modulus in the subcode CELMOD. Although FRAPCON-3.4 includes the coding to return the 

standard error of the model, neither FRAPCON-3.4 nor FRAPTRAN 1.4 uses these errors. 

FRAPTRAN 1.4 includes additional coding to determine Young’s modulus for ZrNb-1. 

CSHEAR

In MATPRO, the expressions used to determine shear modulus in the alpha and beta phases are 

given in Equations 3.6-11 and 3.6-12, respectively. 

 

(3.6-11) 

(3.6-12) 

In the alpha plus beta phase, MATPRO determines G by linearly interpolating between values 

calculated at the alpha to alpha plus beta and the alpha plus beta to beta phase boundaries. 

1 3 2

11 7     c/)cwkf*cdeloxy*cctemp*10x475.510x088.1(modcel

 10x16.1c  7015.5*10x037.1*ctemp 11 8

1  

0.1c2 

10

3  10x6.2c

     25 25

2   10x1/fnckexp10x1/fnckexp1*88.0c

 

3

21

10 7

K

KKT10x168.210x04.4 G   

T10x66.110x49.3G 10 7  
3.21 

The constants K1, K2, and K3 have the same definitions as stated before in CELMOD. However, 

the expression used to model the effect of oxidation for shear modulus is given in Equation 

3.6-13. 

(3.6-13) 

The standard error of the CSHEAR code is 9x109

 Pa. 

The CSHEAR subcode included in the FRAP codes differs slightly from the MATPRO version. 

For instance, the expression for shear modulus in the alpha phase is presented in Equation 3.6-14. 

(3.6-14) 

where 

cshear = shear modulus (Pa) 

ctemp = cladding temperature (K) 

deloxy = input average oxygen concentration excluding oxide layer 

cwkf = input effective cold work (unitless ratio of areas) 

In Equation 3.6-14, c1, c3, and c2 are expressions that account for oxygen content, cold work, and 

fast neutron fluence, respectively. 

Equations 3.6-15 through 3.6-17 describe the constants c1, c2, and c3, respectively. 

(3.6-15) 

(3.6-16) 

(3.6-17) 

For neutron fluences greater than 1x1022, c2 is given by Equation 3.6-18. 

(3.6-18) 

where 

fnck = input effective fast fluence (n/m2

) 

Both FRAPCON-3.4 and FRAPTRAN 1.4 use the same correlations to determine Young’s 

modulus in the subcode CSHEAR. Although FRAPCON-3.4 includes the coding to return the 

standard error of the model, neither FRAPCON-3.4 nor FRAPTRAN 1.4 returns the error. 

FRAPTRAN 1.4 includes additional coding to determine Young’s modulus, shear modulus, and 

Poisson’s ratio as a function of temperature for ZrNb-1. 

Since there is limited data available from shear modulus measurements, a model-to-data 

comparison for Young’s modulus is presented. Figure 3.6-1 presents predicted vs. measured 

Young’s modulus values for data sets used in the development of the CELMOD models as well 

as more recent data. 

  )T10x315.210x07.7(K  11 8

1

1 3 2

10 7     c/)cwkf*cdeloxy*cctemp*10x168.210x04.4(cshear

11 8

1  10x07.7c  10x315.2*ctemp

0.1c2 

10

3  10x867.0c

     25 25

2   10x1/fnckexp10x1/fnckexp1*88.0c
3.22 

The upper and lower bounds of Figure 3.6-1 represent an interval of ± 6.4 GPa. It may be noted 

that the data fit is slightly better below 80 GPa ( ~ < 4 percent) than above 80 GPa 

(~ > 5 percent). This is attributed to the relatively low temperatures that were used to collect 

these data (e.g., a quarter of Bolmaro’s data is collected at < 250°F). Nevertheless the CELMOD 

model, in general, provides a satisfactory fit to the data within the temperature range it will be 

applied for predicting in-reactor performance (> 250°F). 
3.23 

0

4080

120

0 40 80 

Measured Young's Modulus, GPa

Predicted Young's Modulus, GPa

Lower Bound Upper Bound Armstrong and BroPadel and Groff Busby SpasicMehan Northwood Bolmaro and PovoMeasured = Predicted

Figure 3.6-1 Predicted vs. measured values for Young’s modulus for various data sets. 
3.24 

3.6.2 Applicability and Uncertainty 

The CELMOD and CSHEAR correlations are applicable below the melting temperature of the 

cladding (2098K). Above this temperature, the moduli are assumed to be zero (but actually set at 

1x10-10 to avoid dividing by zero). The CELMOD and CSHEAR correlations are similar between 

MATPRO and the FRAPCON-3.4/FRAPTRAN 1.4. The uncertainty of the predictions for 

CELMOD are within 6.4 GPa (the upper and lower bounds), and typically < 5 percent below 

80 GWd/MTU. The uncertainty of the CSHEAR predictions is estimated to be within 9 GPa. 

3.6.3 References 

Bunnell, L.R., et al. 1977. High Temperature Properties of Zircaloy Oxygen Alloys, EPRI NP524, Battelle Pacific Northwest Laboratory, Richland, WA. 

Fisher, E.S., and C.J. Renken. 1964. “Single – Crystal Elastic Moduli and the CHP-BCC 

Transformation in Ti, Zr, and Hf,” Physical Review, 135(2A):A482-494. 

Armstrong, P.E., and H.L. Brown. 1964. “Dynamic Young’s Modulus Measurement above 

1000°C on Some Pure Polycrystalline Metals and Commercial Graphites,” Transactions of the 

Metallurgical Society of AIME 230, August 1964, pp. 962-966. 

Padel, A., and A. Groff. 1976. “Variation du Module de Young du Zirconium b en Function de 

la Temperature,” Journal of Nuclear Materials 59:325-326. 

Shober, F.B., et al. 1957. The Mechanical Properties of Zirconium and Zircaloy-2, BMI-1168, 

Battelle Memorial Institute, Columbus, OH. 

Whitmarsh, C.L. 1962. Review of Zircaloy-2 and Zircaloy-4 Properties Relevant to N.S. 

Savannah Reactor Design, ORNL-3281, Oak Ridge National Laboratory, Oak Ridge, TN. 

Mehan, R.L. 1958. Modulus of Elasticity of Zircaloy-2 Between Room Temperature and 

1,000°F, KAPL-M-RLM-16, General Electric Company, Knolls Atomic Power Laboratory, 

Schenectady, NY. 

Mehan, R.L., and F.W. Wiesinger. 1961. Mechanical Properties of Zricaloy-2, KAPL-2110, pp. 

11-12, General Electric Company, Knolls Atomic Power Laboratory, Schenectady, NY. 

Busby, C.C. 1966. Properties of Zircaloy-4 Tubing, WAPD-TM-585, p. 65, Bettis Atomic 

Power Laboratory, Pittsburgh, PA. 

Northwood, D.O., et al. 1975. “Elastic Constants of Zirconium Alloys,” Journal of Nuclear 

Materials, 55:299-310. 

Spasic, Z., et al. 1968. Conference on the Use of Zirconium Alloys in Nuclear Reactors, 

Marlanske Lanze, Czechoslovakia, USAEC CONF-681086, 1968, pp. 277-284. 

Bolmaro, R.E., and F. Povolo. 1988. “Anelastic behaviour of materials under multiaxial strains,” 

J. Materials Science, 23:371. 
3.25 

3.7 Cladding Axial Growth (CAGROW) 

The subcode CAGROW calculates the factional change in length of Zircaloy tubes due to 

irradiation-induced growth. The change in length of commercial fuel rods due to irradiation 

growth is small. However, growth can be a significant fraction of the clearance between the rod 

and the top and bottom assembly nozzles. Contact with the nozzles can cause rods to bow and 

possibly fail at points where the rods contact each other. In addition, rod growth increases the 

internal void volume that impacts the rod internal pressure calculation. 

The original CAGROW model used in MATPRO was modified in FRAPCON-3.4 to better 

describe the cladding axial growth at higher fluence. The MATPRO version of CAGROW was 

applied to both Zircaloy-2 and -4. The FRAPCON-3.4 version includes correlations for 

Zircaloy-2 and -4, M5, and ZIRLO. These correlations may be applied to either PWR or boiling 

water reactor (BWR) conditions. The CAGROW subroutine is not used in FRAPTRAN 1.4 

because significant irradiation axial growth is not expected under the short time scales considered 

by FRAPTRAN 1.4. 

3.7.1 Model Development and Comparisons 

The MATPRO correlation used to describe cladding axial growth is presented in Equation 3.7-1. 

This equation has been developed to model the irradiation growth of Zircaloy tubes at 

temperature between 40 and 360°C (the normal range of cladding temperatures in LWRs). 

(3.7-1) 

where 

L/L = fractional change in length due to growth. 

A = 1.407x10-16 (n/m2

)

1/2.

T = cladding temperature (K). 

 = fast neutron flux (n/m2

s) (E > 1.0 MeV). 

t = time (s). 

fz = texture factor for the tubing. The variable fz is the effective fraction of cells 

aligned with their <0001> axis parallel to the tubing axis, as determined by Xray diffraction analysis. A value of fz = 0.05 is typical. 

CW = cold work (fraction of cross-sectional area reduction). 

Axial growth for temperature below 40°C is approximated by using T = 40°C in Equation 3.7-1 

and growth above 360°C is approximated by using T = 360°C. 

The MATPRO subroutine CAGROW was based on low-burnup, low-fluence data found in 

Harbottle (1970), Kreyns (1960), and Daniel (1971, 1972). Consequently, the CAGROW model 

was updated in FRAPCON-3.4 to better model axial growth under greater fluences ( > 1x1025

n/m2

). A model proposed by Franklin (1982) that was based on high fluence PWR data was 

adopted for FRAPCON-3.4. This model is presented in Equation 3.7-2 and relates the axial 

growth, ax (m/m), of stress relief annealed (SRA) Zircaloy-4 in a PWR to fluence (n/cm2

). 

(3.7-2) 

        CW02.01f1tT/8.240expAL/L  z

2/1  

845.021  *10x18.2ax  
3.26 

where 

ax = axial growth increment (m/m) 

 = fast neutron fluence (n/cm2

) (E > 1.0 MeV) 

The expression in Equation 3.7-2 is used to determine the axial growth at the beginning and the 

end of a time increment. The difference is used to determine the axial growth for the time 

increment, which is the value returned by the CAGROW subroutine. The expression in Equation 

3.7-2 can be used to describe the axial growth of fully recrystallized annealed (RXA) Zircaloy-2 

cladding in a BWR by multiplying by a factor of 0.5. 

In addition to updating the cladding axial growth model for SRA Zircaloy-4 and RXA Zircaloy-2, 

new correlations have been added to describe the cladding axial growth of M5 and ZIRLO. The 

cladding axial growth correlations for M5 and ZIRLO are presented in Equations 3.7-3 and 3.7-4, 

respectively. These equations relate the cladding axial growth, ax (m/m), to fluence (n/cm2

). The 

variables have the same definitions as those described in Equation 3.7-2 and the time difference is 

used to determine the axial growth during the time increment, which is the value returned by the 

CAGROW subroutine. 

For M5 cladding: 

(3.7-3) 

For ZIRLO cladding: 

 (3.7-4) 

Model-to-data comparisons are presented in Figures 3.7-5 through 3.7-8 for SRA Zircaloy-4, 

RXA Zircaloy-2, M5, and ZIRLO, respectively. Based on these figures, it appears that the 

cladding axial growth of these alloys can be predicted well with the CAGROW subroutine. 

81787.021 10013.7ax  

98239.025 107893.9ax  
3.27 

0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

0 2E+21 4E+21 6E+21 8E+21 1E+22 1.2E+22

Fast Neutron Fluence, n/cm2-s

L/L, %

Newman Franklin Gilbon Zr-4 PWR Model

Figure 3.7-5 Model-to-data comparison for SRA Zircaloy-4 (PWR). The standard error of 

L/L (%) is 0.11. 

0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

0 2E+21 4E+21 6E+21 8E+21 1E+22 1.2E+22

Fast Neutron Fluence, n/cm2-s

L/L, %

FRAPCON Data Gilbon Zr-2 BWR Model

Figure 3.7-6 Model-to-data comparison for RXA Zircaloy-2 (BWR). The standard error 

of L/L (%) is 0.07. 
3.28 

0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

0 2E+21 4E+21 6E+21 8E+21 1E+22 1.2E+22

Fast neutron fluence, n/cm2, E>1MeV

Fuel rod growth, %

M5 Data PWR M5 Model

Figure 3.7-7 Model-to-data comparison for M5. M5 data was collected by Gilbon et al. 

(2000). The standard error of L/L (%) is 0.05. 

0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

0 2E+21 4E+21 6E+21 8E+21 1E+22 1.2E+22

Fast neutron fluence, n/cm2, E>1MeV

Fuel rod growth, %

Vandellos ZIRLO Data North Anna ZIRLO Data

BR-3 ZIRLO Data PWR ZIRLO Model

Figure 3.7-8 Model-to-data comparison for ZIRLO. Vandellos data was collected by Irisa 

(2000) while North Anna and BR-3 data was collected by Sabol (1994). The 

standard error of L/L (%) is 0.05. 
3.29 

3.7.2 Applicability and Uncertainty 

The CAGROW model in FRAPCON-3.4 is applicable over the temperature range from 700 to 

900K for local burnups up to 65 GWd/MTU (or 12x1025 n/m2 fluence, E > 1MeV). The absolute 

uncertainties for the growth strains predicted for stress-relief annealed Zr-4, M5, and ZIRLO are 

0.11, 0.05, and 0.05, respectively. These three alloys are typically used in PWRs, whereas 

recrystallized Zr-2 is used in BWRs. The absolute uncertainty in growth strains for Zr-2 is 0.07. 

3.7.3 References 

Harbottle, J.E. 1970. “The Temperature and Neutron Dose Dependence of Irradiaton Growth in 

Zircaloy-2,” Irradiation Effect on Structural Alloys for Nuclear Reactor Applications, ASTMSTP-485, American Society of Testing and Materials, West Conshohocken, PA, pp. 287-299. 

Kreyns, P.H., quoted by E. Duncombe et al. 1996. Comparison with Experiment of Calculated 

Changes and Failure Analysis of Irradiated Bulk Oxide Fuel Test Rods Using the CYGRO01 

Computer Program, WAPD-TM-583, Bettis Atomic Power Laboratory, Pittsburgh, PA. 

Daniel, R.C. 1971. In-Pile Dimensional Changes of Zircaloy-4 Tubing Having Low Hoop 

Stresses (LWBR Development Program), WAPD-TM-973, Bettis Atomic Power Laboratory, 

Pittsburgh, PA. 

Daniel, R.C. 1972. “In-Pile Dimensional Changes of Zircaloy-4 Tubing Having Stresses (Light 

Water Breeder Reactor Development Program),” Nuclear Technology, 14:171-186. 

Franklin, D.G. 1982. “Zircaloy Cladding Deformation during Power Reactor Irradiation.” In 

Proceedings of the Fifth International Symposium on Zirconium in the Nuclear Industry, ASTMSTP-754, American Society of Testing and Materials, West Conshohocken, PA, pp. 235-267. 

Newman, L.W. 1986. The Hot Cell Examination of Oconee 1 Fuel Rods After Five Cycles of 

Irradiation, DOE/ET/34212-50, Babcock and Wilcox Company, Lynchburg, VA. 

Gilbon, D., A. Soniak, S. Doriot, and J.P. Mardon. 2000. “Irradiation Creep and Growth 

Behavior, and Microstructural Evolution of Advanced Zr-Base Alloys,” Zirconium in the Nuclear 

Industry: Twelfth International Symposium, ASTM-STP-1354, GP Sabol and GD Moan eds., 

American Society of Testing and Materials, West Conshohocken, PA, pp. 51-73. 

Irisa, Y., Murai, K., Goto, K., Suzuki, S., Sabate, R., Alonso, J.M. 2000. “Segmented Fuel Rod 

Irradiation Programme on Advanced Materials for High Burnup,” in Proceedings of the 2000 

Light Water Fuel Performance Meeting, Park City, UT, April 10-13, 2000. 

Sabol, G.P., R.J. Comstock, R.A. Weiner, P. Larouere, and R.N. Stanutz. 1994. “In-Reactor 

Corrosion Performance of ZIRLO and Zircaloy-4,” Zirconium in the Nuclear Industry: Tenth 

International Symposium, ASTM-STP-1425, A.M. Garde and E.R. Bradley, eds. American 

Society for Testing and Materials, Philadelphia, PA, pp. 724-744. 
3.30 

3.8 Creep Rate (CREPR) 

The subcode CREPR is used by FRAPCON-3.4 to determine the cladding creep rate, accounting 

for both thermal and irradiation creep. The CREPR subcode was not used in MATPRO. Also, 

since creep is a time-dependent deformation process, CREPR is not used in FRAPTRAN 1.4, 

because creep occurs over much longer time intervals than those considered in FRAPTRAN 1.4 

of a minute or less. The FRAPTRAN calculation of cladding ballooning as a result of a LOCA is 

handled with a separate subroutine named BALON2 

3.8.1 Model Development and Comparisons 

The creep model used by FRAPCON-3.4 is based on the model given by Limbäck and Andersson 

(1996). This model uses a thermal creep model described by Matsuo (1987) and an empirical 

irradiation creep rate with tuned model parameters that were fit to data by Franklin et al. (1983). 

The Limbäck model was further modified by PNNL to use effective stress rather than hoop stress 

as an input so that the difference in creep behavior during tensile and compressive creep would be 

modeled correctly. Several of the fitting coefficients were consequently changed to 

accommodate this modification. In addition, a temperature-dependent term was added to the 

formula for irradiation creep strain rate. This model has different parameters for SRA and RXA 

cladding types, and provides much more reasonable creep strains in the LWR range of 

temperatures and cladding hoop stresses than the previous model used in FRAPCON-3.3. 

The steady-state thermal and irradiation creep rates are given by Equations 3.8-1 and 3.8-2, 

respectively. 

(3.8-1) 

(3.8-2) 

where 

th  = thermal strain rate (in./in./hr) 

irr  = irradiation strain rate (in./in./hr) 

T = temperature (K) 

eff = effective stress (MPa) 

R = universal gas constant (0.008314 kJ/mol-K) 

 = fast neutron flux (n/m2

-s) 

The variables A, E, ai, n, -Q, Co, C1, C2, and the function f(T) have different values depending on 

temperature, flux, and cladding type (SRA or RXA). The values for these variables under 

different conditions are presented in Table 3.8-1. 





 



  



 



 



RT

Q exp Ea sinh

T

E A

n

effi

th 

C )T(f 1 C2

eff

C

oirr  
3.31 

Table 3.8-1 Parameters for FRAPCON-3.4 creep equation for SRA and RXA cladding. 

Parameter Units SRA Cladding RXA Cladding 

A* K/MPa/hr 1.08x109* 5.47x108* 

E MPa 1.149x105 -59.9*T 

ai MPa-1 650{1-0.56[1-exp(-1.4x10-27*1.3)]} 

=fast neutron fluence (n/cm2

) 

n unitless 2.0 3.5 

Q kJ/mole 201 

Co 

(n/m2-s)-Cl 

MPa-C2 4.0985x10-24 1.87473x10-24 

C1 unitless 0.85 

C2 unitless 1.0 

f(T)* unitless 

T < 570 0.7283 

570< T< 625 -7.0237+0.0136T 

T > 625 1.4763 

0.7994 

-3.18562+0.00699132T 

1.1840 

*Parameters changed from original Limbäck equation. 

The thermal and irradiation creep rates may be added together as shown in Equation 3.8-3 and 

used to calculate the saturated primary hoop strain, which is presented in Equation 3.8-4. 

(3.8-3) 

(3.8-4) 

The total thermal strain can then be calculated as a function of time, t (hours), as shown in 

Equation 3.8-5. However, in FRAPCON-3.4 the strain rate is used, which is obtained by taking 

the derivative of Equation 3.8-5. This derivative is presented in Equation 3.8-6, which relates the 

total thermal strain rate to the saturated primary hoop strain, the combined thermal and irradiation 

creep rates, and time, t, in hours. 

(3.8-5) 

(3.8-6) 

The effective stress in the cladding is found using the principal stresses at the mid-wall radius 

using the thick wall formula. The principle stresses can be determined with Equations 3.8-7 

through 3.8-9. 

(3.8-7) 

(3.8-8) 

irrthirrth       

   05.2

irrth

109.0

irrth

s

p 0216.0 35500tanh2 

    



  52exp1  tt irrth irrth

s

pH    

  2/1 irrth irrth

2/1

irrth

s

p

H 52exp t

t2

52

 

  



   

 

2

i

2

o

2

io

2

o

2

2 i

oo

2

ii

r rr

r

PPrr rPrP



 



 

2

i

2

o

2

io

2

o

2

2 i

oo

2

ii

t rr

r

PPrr rPrP



 


3.32 

(3.8-9) 

where 

Pi = inner pressure 

Po = outer pressure 

ri = inner radius 

ro = outer radius 

r = radius within tube 

r = radial stress 

t = tangential stress 

l = longitudinal stress 

The effective stress can then be calculated by Equation 3.8-10. 

(3.8-10) 

The subcode CREPR can be used for newer alloys like M5 and ZIRLO. It has been found that 

the Zircaloy RXA model adequately describes the creep behavior of M5 (Gilbon et al., 2000; 

Soniak et al., 2002). The Zircaloy SRA model is used for ZIRLO with a reduction factor of 0.8. 

The reduction factor is the result of studies that have shown that ZIRLO exhibits about 80 percent 

of the Zircaloy-4 creepdown (Sabol et al., 1994). 

Model-to-data comparisons are presented in Figures 3.8-1 and 3.8-2 for SRA and RXA Zircaloy 

cladding. Creep data was obtained from irradiated SRA and RXA tube from Franklin (1983), 

Soniak et al. (2002), and Gilbon et al. (2000). The relative standard error for the predictions is 

14.5 percent for SRA and 21.6 percent for RXA cladding. 

2

i

2

o

2

oo

2

ii

l

rr

rPrP



 

     2

lr

2

rt

2

eff 5.0 tl 
3.33 

0

0.005

0.01

0.015

0.02

0.025

0 0.005 0.01 0.015 0.02 0.025

Measured hoop strain, in/in

Predicted hoop strain, in/in

Franklin Soniak 350°C Soniak 320°C Gilbon Predicted=Measured

Figure 3.8-1 Predicted vs. measured hoop strain for SRA Zircaloy model. 

0

0.002

0.004

0.006

0.008

0 0.002 0.004 0.006 0.008

Measured hoop strain, in/in

Predicted hoop strain, in/in

Franklin Soniak 350°C Soniak 320°C Gilbon Predicted=Measured

Figure 3.8-2 Predicted vs. measured hoop strain for RXA Zircaloy model. 

3.8.2 Applicability and Uncertainty 

The subroutine CREPR is applicable over the range of conditions for which the data was 

collected. The database spans a temperature range of 570 to 625K, an effective stress range of 40 

to 130 MPa, and a fast neutron flux range of 1x1017 to 2x1018 n/m2

-s. This covers the nominal 
3.34 

operating conditions for LWRs. The relative standard error for the predictions is 14.5 percent for 

SRA and 21.6 percent for RXA cladding. The uncertainty for predicting cladding creep in a fuel 

rod will most likely be greater than these values because there may be greater uncertainty in the 

estimate of fast neutron flux and cladding temperature than those from the creep data used to 

determine the creep model uncertainties. 

3.8.3 References 

Limbäck, M., and T. Andersson. 1996. “A Model for Analysis of the Effect of Final Annealing 

on the In- and Out-of-Reactor Creep Behavior of Zircaloy Cladding,” Zirconium in the Nuclear 

Industry: Eleventh International Symposium, ASTM STP 1295, E.R. Bradley and G.P. Sabol, 

eds., American Society for Testing and Materials, West Conshohocken, PA, pp. 448-468. 

Matsuo, Y. 1987. “Thermal Creep of Zircaloy-4 Cladding under Internal Pressure,” Journal of 

Nuclear Science and Technology, 24(2):111-119. 

Franklin, D.G., G.E. Lucas, A.L. Bement. 1983. “Creep of Zirconium Alloys in Nuclear 

Reactors,” ASTM STP 815, American Society for Testing and Materials, West Conshohocken, 

PA. 

Gilbon, D., A. Soniak, S. Doriot, and J.P. Mardon. 2000. “Irradiation Creep and Growth 

Behavior, and Microstructural Evolution of Advanced Zr-Base Alloys,” Zirconium in the Nuclear 

Industry: Twelfth International Symposium, ASTM-STP-1354, G.P. Sabol and G.D. Moan, eds., 

American Society of Testing and Materials, West Conshohocken, PA, pp. 51-73. 

Soniak, A., N. L’Hullier, J-P. Mardon, V. Rebeyrolle, P. Boffioux, C. Bernaudat. 2002. 

“Irradiation Creep Behavior of Zr-Base Alloys,” Zirconium in the Nuclear Industry: Thirteenth 

International Symposium, ASTM STP 1423, G.D. Moan and P. Rudling, eds., ASTM 

International, West Conshohocken, PA, pp.837-862. 

Sabol, G.P., Comstock, R.J., Weiner, R.A., Larouere, P., Stanutz, R.N. 1994. “In-Reactor 

Corrosion Performance of ZIRLOTM and Zircaloy-4,” Zirconium in the Nuclear Industry: Tenth 

International Symposium, ASTM STP 1245, A.M. Garde and E.R. Bradley, eds. American 

Society for Testing and Materials, Philadelphia, PA, pp. 724-744. 

3.9 Cladding Meyer Hardness (CMHARD) 

The subroutine CMHARD calculates Meyer hardness as a function of cladding temperature. 

Hardness is one of the parameters required for calculating fuel-to-cladding contact conductance. 

As the contact pressure between the two surfaces increases, the points of contact enlarge due to 

localized plastic deformation and the solid-to-solid thermal conductance is improved. The Meyer 

hardness is used by Ross and Stoute (1962) in their heat transfer correlation as an indication of 

the hardness of resistance to deformation of the softer (Zircaloy) material. 

The same CMHARD subroutine is used in the MATPRO, FRAPCON-3.4, and FRAPTRAN 1.4 

codes. However, FRAPTRAN 1.4 includes additional coding that ensures that the minimum 

hardness returned is 1.94x108

 N/m2 (the highest temperature data point) and includes provisions 

for ZrNb-1. 
3.35 

3.9.1 Model Development and Comparisons 

In MATPRO, the Meyer hardness number is a measure of indentation hardness and is defined in 

conjunction with Meyer’s law, which is presented in Equation 3.9-1. 

(3.9-1) 

where 

L = load 

d = the diameter of impression at the surface of a specimen in a static ball test 

n = the Meyer work hardening coefficient 

a = a material constant 

The Meyer hardness number (MH) is defined as 4L/d2

. Other hardness numbers are available 

(Brinell, Rockwell, etc.), and conversion from one to another is possible. However, the routine 

CMHARD was created to provide information required by the Ross and Stoute gap conductance 

model that includes a dependence on Meyer hardness. 

Meyer hardness numbers for temperatures from 298 to 877K were taken from Peggs and Godin 

(1975). A regression analysis of the reciprocal of the Meyer hardness values versus the log of 

temperature was used to obtain the analytical expression used in CMHARD. The correlation 

used is given by Equation 3.9-2. 

(3.9-2) 

where 

MH = Meyer hardness (N/m2

) 

T = temperature (K) 

The Meyer hardness decreased rapidly with increasing temperature, beginning at 2x109

 MPa at 

room temperature and decreasing to 2x108

 MPa at 875K. The hardness is presumed to continue 

its rapid rate of decrease at temperatures above 875K. The minimum Meyer hardness number of 

Zircaloy cladding is 1.0x105

 N/m2. 

A comparison between the CMHARD correlation and the data from which it was derived is 

presented in Figure 3.9-1. 

n  adL

     1 2 5 8 10x5621.2T10x3502.4T10x6394.2T10x6034.2expMH       
3.36 

0

200

400

600

800

1000

1200

1400

1600

1800

2000

200 400 600 800 1000 1200 1400 1600 1800 2000 2200

Temperature (K)

Meyer Hardness (MPa)

Figure 3.9-1 CMHARD correlation as a function of temperature with data from Peggs and 

Godin. 

3.9.2 Applicability and Uncertainty 

The subroutine CMHARD is applicable for temperatures up to 875K. An estimate of the 

uncertainty has not been established for this correlation due to the limited data. 

3.9.3 References 

Ross, A.M., and R.L. Stoute. 1962. Heat Transfer Coefficient Between UO2 and Zircaloy-2, 

AECL-1552, Atomic Energy of Canada, Ltd., Chalk River, Ontario. 

Peggs, I.D., and D.P. Godin. 1975. “The Yield Strength-Hot Hardness Relationship of Zircaloy4,” Journal of Nuclear Materials, 57:246-248. 
4.1 

4.0 Gas Material Properties 

This section describes material property correlations for gap gases. The subroutine used to 

determine the material correlation in either the FRAPCON-3/FRAPTRAN codes is given in 

parentheses in the heading of each subsection. 

4.1 Gas Conductivity (GTHCON) 

The subroutine GTHCON calculates the gas thermal conductivity as a function of temperature 

and gas fraction for seven gases: helium (He), argon (Ar), krypton (Kr), xenon (Xe), hydrogen 

(H), nitrogen (N), and water vapor (steam). The MATPRO, FRAPCON-3.4, and FRAPTRAN 

1.4 codes use similar correlations to determine the gas thermal conductivity. However, 

FRAPCON-3.4 and FRAPTRAN 1.4 use updated fitting parameters to better estimate gas thermal 

conductivity at higher temperatures. 

4.1.1 Model Development and Comparisons 

The heat conductance of gas-filled gaps or pores is dependent on the thermal conductivity of the 

gas mixture when the dimensions of the gas-filled regions are large compared to the mean 

distance between gas molecule collisions (mean free path of the gas molecules). This section 

presents data and correlations for the thermal conductivities of the gases of interest in fuel rod 

analysis. 

The correlations used for the gases of interest are all of the form presented in Equation 4.1-1. 

(4.1-1) 

where 

K = thermal conductivity (W/m-K) 

T = gas temperature (K) 

The constants A and B are fitting parameters. The values for A and B used in MATPRO for each 

gas are given in Table 4.1-1 and the values used in FRAPCON-3.4 and FRAPTRAN 1.4 (FRAP) 

are given in Table 4.1-2. 

Table 4.1-1 Constants used in gas thermal conductivity correlations in MATPRO. 

Gas A B 

He 2.639x10-3 0.7085 

Ar 2.986x10-4 0.7224 

Kr 8.247x10-5 0.8363 

Xe 4.351x10-5 0.8616 

H2 1.097x10-3 0.8785 

N2 5.314x10-4 0.6898 

B K  AT
4.2 

Table 4.1-2 Constants used in gas thermal conductivity correlations in FRAP. 

Gas A B 

He 2.531x10-3 0.7146 

Ar 4.092x10-4 0.6748 

Kr 1.966x10-4 0.7006 

Xe 9.825x10-5 0.7334 

H2 1.349x10-3 0.8408 

N2 2.984x10-4 0.7799 

The MATPRO steam correlation is also used in FRAPCON-3.4/FRAPTRAN 1.4 codes. For 

temperatures less than or equal to 973.15K, the correlation presented in Equation 4.1-2 is used. 

For temperatures greater than 973.15K, the correlation presented in Equation 4.1-3 is used to 

determine the thermal conductivity of steam. Temperature, T, is in kelvin and pressure, P, is in 

pascals. 

(4.1-2) 

(4.1-3) 

The thermal conductivity of gas mixtures is calculated with Equation 4.1-4. 

(4.1-4) 

where ij is given in Equation 4.1-5 

(4.1-5) 

and ij is given in Equation 4.1-6 

(4.1-6) 

and 

ij = Kronecker delta = 1 for i = j, 0 otherwise (unitless) 

n = number of components in mixture (unitless) 

Mi = molecular weight of component i (kg) 

 

   

  ]15.273T*10x51.415.273T*10x08.1

15.273T*10x87.510x6.17[15.273T/T/P*009.1

T*10x005.6T*10x424.910x8516.2*T/Pk

7 2 11 3

22 2.4 4 5

8 10 214

  

   

   

 

 

  





45.16 5 9 3.1)T/P*10x1668.2(*10x45.9T*10x44.4k     

  

 

































n

i

n

1j

i jijij

ii

mix

x1x

xk k

  

  

















  2

ji

iji j

ijij

MM

M142.0MMM

41.21

2/1

j

2/3 i

2 4/1

j

i

2/1

j

i

ij

M

M 12

M

M

k

k 1






































































4.3 

xi = mole fraction of component i (unitless) 

ki = thermal conductivity of the component i (W/m-K) 

Model-to-data comparisons for the gas conductivity models are presented in Figures 4.1-1 

through 4.1-7. 

0.00

0.10

0.20

0.30

0.40

0.50

0.60

0.70

0.80

0 500 1000 1500 2000 2500 3000

Temperature, K

Thermal Conductivity, W/m-K

MATPRO FRAP

Faubert and Springer (1973) Jain and Saxena (1975)

Jody et al. (1977)

Figure 4.1-1 Model-to-data comparison for He thermal conductivity. 
4.4 

0.00

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

0.10

0 500 1000 1500 2000 2500 3000

Temperature, K

Thermal Conductivity, W/m-K

MATPRO FRAP

Timrot and Umanski (1966) Saxena and Saxena (1968)

Faubert and Springer (1972) Springer and Wingeier (1973)

Stefanov et al. (1976)

Figure 4.1-2 Model-to-data comparison for Ar thermal conductivity. 
4.5 

0.00

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0 500 1000 1500 2000 2500 3000

Temperature, K

Thermal Conductivity, W/m-K

MATPRO FRAP

Saxena and Saxena (1969) Stefanov et al. (1976)

Vargaftik and Yakush (1971)

Figure 4.1-3 Model-to-data comparison for Kr thermal conductivity. 
4.6 

0.00

0.01

0.01

0.02

0.02

0.03

0.03

0.04

0.04

0.05

0 500 1000 1500 2000 2500 3000

Thermal Conductivity, W/m-K

Temperature, K

MATPRO FRAP

Stefanov (1976) Springer and Wingeier (1973)

Saxena and Saxena (1969) Vargaftik and Yakush (1971)

Figure 4.1-4 Model-to-data comparison for Xe thermal conductivity. 

0.00

0.20

0.40

0.60

0.80

1.00

1.20

0 500 1000 1500 2000 2500 3000

Temperature, K

Thermal Conductivity, W/m-K

MATPRO FRAP Saxena and Saxena (1970)

Figure 4.1-5 Model-to-data comparison for H2 thermal conductivity. 
4.7 

0.00

0.02

0.04

0.06

0.08

0.10

0.12

0.14

0.16

0 500 1000 1500 2000 2500 3000

Temperature, K

Thermal Conductivity, W/m-K

MATPRO FRAP

Vargaftik and Zimina (1964) Faubert and Springer (1972)

Chen and Saxena (1973)

Figure 4.1-6 Model-to-data comparison for N2 thermal conductivity. 

0

0.02

0.04

0.06

0.08

0.1

0.12

0.14

0.16

600 700 800 900 1000

Temperature (K)

Thermal Conductivity (W/m-K)

MATPRO/FRAP

Data

Figure 4.1-7 Model-to-data comparison for steam thermal conductivity. 

4.1.2 Applicability and Uncertainty 

The gas thermal conductivity correlations for used in GTHCON are applicable over the 

temperature range for which data was collected. The standard error for each correlation is given 

in Table 4.1-3. 
4.8 

Table 4.1-3 Standard deviation of gas thermal conductivity correlations. 

Gas Standard Deviation 

He 8.99x10-3

Ar 9.66x10-4

Kr 8.86x10-4

Xe 5.34x10-4

H2 1.67x10-2

N2 1.97x10-4

Steam 1.69x10-3

4.1.3 References 

Faubert, F.M., and G.S. Springer. 1973. “Measurement of the Thermal Conductivity of Helium 

up to 2100°K by the Column Method,” Journal of Chemical Physics, 58(10):4080-4083. 

Jain, P.C., and S.C. Saxena. 1975. “Transport Properties of Helium in the Temperature Range 

400-2300K,” Chemical Physics Letters, 36(4):489-491. 

Jody, B.J., S.C. Saxena, V.P.S. Nain, and R.A. Aziz. 1977. “Thermal Conductivity of Helium: 

A Problem for the Repulsive Wall of the Interatomic Potential,” Chemical Physics, 22:53-58. 

Timrot, D.L., and A.S. Umanskii. 1966. “Thermal Conductivity of Hydrogen and Argon,” High 

Temperature, 3:285-287. 

Saxena, V.K., and S.C. Saxena. 1968b. “Measurement of the Thermal Conductivity of Argon 

Using Hot-Wire Type Thermal Diffusion Columns,” Chemical Physics Letters, 2(1):44-46. 

Faubert, F.M., and G.S. Springer. 1972. “Measurement of the Thermal Conductivity of Argon, 

Krypton, and Nitrogen in the Range 800-2000°K,” Journal of Chemical Physics, 57(6):2333-

2340. 

Springer, G.S., and E.W. Wingeier. 1973. “Thermal Conductivity of Neon, Argon, and Xenon at 

High Temperatures,” Journal of Chemical Physics, 59(5):2747-2750. 

Stefanov, B., L. Zarkova, and D. Oliver. 1976. “Measurement of the Thermal-Conductivity 

Coefficient of Gases and Vapors up to 2500°K,” High Temperature, 14(1):48-57. 

Saxena, V.K., and S.C. Saxena. 1969. “Thermal Conductivity of Krypton and Xenon in the 

Temperature Range 350-1500K,” Journal of Chemical Physics, 51(8):3361-3368. 

Vargaftik, N.B., and L.V. Yakush. 1971. “Measurement of the Thermal Conductivities of Neon, 

Krypton, and Xenon over a Wide Range of Temperatures,” Journal of Engineering Physics, 

21(3):1156-1161. 

Saxena, S.C., and V.K. Saxena. 1970. “Thermal Conductivity Data for Hydrogen and Deuterium 

in the Range 100-1100°C,” Journal of Physics A, 3(3):309-320. 

Vargaftik, N.B., and N.Kh. Zimina. 1964. “Thermal Conductivity of Nitrogen at High 

Temperatures,” High Temperature, 2(6):782-790. 
4.9 

Chen, S.H.P., and S.C. Saxena. 1973. “Experimental Determination of Thermal Conductivity of 

Nitrogen in the Temperature Range 100-2200°C,” High Temperature Science, 5(3):206-233. 
Appendix A

Subroutine Source Codes 
A.1 

Appendix A Subroutine Source Codes

A.1 Uranium Dioxide/Mixed Oxide Properties 

A.1.1 Fuel Melting Temperature (PHYPRP) 

FRAPCON-3.4: PHYPRP

*deck phyprp 

 subroutine phyprp 

c 

c phyprp returns uo2, (u,pu)o2, and zircaloy melting points 

c and heats of fusion, and zirconium and zircaloy alpha to beta 

c transition temperatures. 

 common / phypro / ftmelt,fhefus,ctmelt,chefus,ctranb, 

 + ctrane,ctranz,fdelta,bu,comp,deloxy 

c ftmelt = output uo2 or mixed oxide fuel melting points (k) 

c fhefus = output uo2 or mixed oxide fuel heat of fusion (j/kg) 

c ctmelt = output zr clad melting point (k) 

c chefus = output zr clad heat of fusion (j/kg). 

c ctranb = output start of zr-4 alpha-beta transus temperature (k) 

c ctrane = output end of zr-4 alpha-beta transus temperature (k) 

c ctranz = output zr isothermal alpha-beta transus temperature (k) 

c fdelta = output liquid-solid coexistence temperature range (k) 

c bu = input burnup (mw-s/kg-u) 

c comp = input puo2 content (wt%) 

c deloxy = input oxygen concentration - oxygen concentration 

c of as received cladding (kg oxygen/kg zircaloy) 

c the uo2 and zr cladding melting points and heats of fusion are 

c taken from the following references. uo2 fuel melting point of 

c 3113k from h.c.brassfield et al gemp-482. uo2 heat of fusion 

c of 17.7 kcal/mole from l.leibowitz et al, j.nuc.mat. 39 p 115 

c (1971). cladding melting point of 2098 k from m.l. picklesimer 

c private communication. cladding heat of fusion of 4.9 kcal/mol 

c from brassfield et al, gemp-482. 

c beginning and end of alpha-beta transus are from data in 

c figure iii.33 of anl-76-49 

c isothermal zirconium alpha-beta transus temperature is 1135 k 

c taken from b.lustman & f.kerze "the metallurgy of zirconium" 

c mcgraw-hill book co., new york, 1955 

c mixed oxide melting point was obtained from lyon et al, j. nuc. 

c mat., 22 (1967) p 332 

c phypro was coded by v.f. baston in may 1974 

c modified by c.s. olsen in feb. 1975 

c modified by b.w. burnham in nov. 1977 

c modified by d.l. hagrman in june 1979 

c burnup dependance for fuel melting temperature changed from 3.2K/GWd/tHM 

c to 0.5 K/GWd/tHM by K.J. Geelhood in June 2004 based on 

c S. G. Popov, et al., "Thermophysical Properties of MOX and UO2 Fuels 

c including the Effects of Irradiation," ORNL/TM-2000/351. 

 include 'lacmdl.h' 

 data on / 1 /, 

 + off / 2 /, 

 + locidx / 2 / 

 sldus(c) = 2840.0-5.41395*c+7.468390e-3*c*c 

 liqdus(c) = 2840.0-3.21860*c-1.448518e-2*c*c 

 fbu = bu/86.4 

 if (comp.gt.0.0) go to 100 

 ftmelt = 3113.15-5.0*fbu/10000.0 
A.2 

 fdelta = 1.0e-10 

 go to 110 

100 c1 = comp 

 ftmelt = sldus(c1)+273.15-5.0*fbu/10000.0 

 fdelta = liqdus(c1)-sldus(c1)-5.0*fbu/10000.0 

110 fhefus = 27.4e+4 

 ctmelt = 2098.15 

 chefus = 22.5e+04 

 wfox = deloxy+0.0012 

 ctranb = 1094.+wfox*(-1.289e+03+wfox*7.914e+05) 

 if (wfox.lt.0.025) go to 120 

 ctranb = 1556.4+3.8281e+04*(wfox-0.025) 

120 ctrane = 392.46*((100.*deloxy+0.1242807)**2+3.1417) 

 if (deloxy.lt.4.7308937e-03) go to 130 

 ctrane = (100.*deloxy+0.12)*491.157+1081.7413 

130 continue 

 ctranz = 1135.15 

 return 

 end 

FRAPTRAN 1.4: PHYPRP

*deck phyprp 

c 

 subroutine phyprp (icm) 

c 

 implicit real (a-h,o-z) 

c 

 real liqdus 

c 

c This is a MATPRO-11, Rev. 2 routine modified by PNNL for use in FRAPT 

c 

c phyprp returns uo2, (u,pu)o2, and zircaloy melting points 

c and heats of fusion, and zirconium and zircaloy alpha to beta 

c transition temperatures. 

c 

 include 'phypro.h' 

c 

c ftmelt = output uo2 or mixed oxide fuel melting points (k) 

c fhefus = output uo2 or mixed oxide fuel heat of fusion (j/kg) 

c ctmelt = output zr clad melting point (k) 

c chefus = output zr clad heat of fusion (j/kg). 

c ctranb = output start of zr-4 alpha-beta transus temperature (k) 

c ctrane = output end of zr-4 alpha-beta transus temperature (k) 

c ctranz = output zr isothermal alpha-beta transus temperature (k) 

c fdelta = output liquid-solid coexistence temperature range (k) 

c 

c bumtp = input burnup (mw-s/kg-u) 

c compmt = input puo2 content (wt%) 

c deloxy = input oxygen concentration - oxygen concentration 

c of as received cladding (kg oxygen/kg zircaloy) 

c 

c the uo2 and zr cladding melting points and heats of fusion are 

c taken from the following references. uo2 fuel melting point of 

c 3113k from h.c.brassfield et al gemp-482. uo2 heat of fusion 

c of 17.7 kcal/mole from l.leibowitz et al, j.nuc.mat. 39 p 115 

c (1971). cladding melting point of 2098 k from m.l. picklesimer 

c private communication. cladding heat of fusion of 4.9 kcal/mol 

c from brassfield et al, gemp-482. 

c 

c beginning and end of alpha-beta transus are from data in 

c figure iii.33 of anl-76-49 

c isothermal zirconium alpha-beta transus temperature is 1135 k 
A.3 

c taken from b.lustman & f.kerze "the metallurgy of zirconium" 

c mcgraw-hill book co., new york, 1955 

c 

c mixed oxide melting point was obtained from lyon et al, j. nuc. 

c mat., 22 (1967) p 332 

c 

c phypro was coded by v.f. baston in may 1974 

c modified by c.s. olsen in feb. 1975 

c modified by b.w. burnham in nov. 1977 

c modified by d.l. hagrman in june 1979 

c Modified by PNNL, January 1997, to clean up coding and delete 

c sensitivity uncertainty analysis coding 

c burnup dependance for fuel melting temperature changed from 

c 3.2 K/GWd/tHM to 0.5 K/GWd/tHM by K.J.Geelhood in June 2004 based on 

c S.G. Popov, et al., "Thermophysical Properties of MOX and UO2 Fuels 

c including the Effects of Irradiation," ORNL/TM-200/351 

c 

c 

 sldus(c) = 2840.0d0 - 5.41395d0*c + 7.468390d-3*c*c 

 liqdus(c) = 2840.0d0 - 3.21860d0*c - 1.448518d-2*c*c 

 fbu = bumtp/86.4d0 

 if (compmt .gt. 0.0d0) go to 10 

c 

 ftmelt = 3113.15d0 - 5.0d0*fbu/10000.0d0 

 fdelta = 1.0d-10 

 go to 20 

c 

 10 continue 

 c1 = compmt 

 ftmelt = sldus(c1) + 273.15d0 - 5.0d0*fbu/10000.0d0 

 fdelta = liqdus(c1) - sldus(c1)- 5.0d0*fbu/10000.0d0 

c 

 20 continue 

 fhefus = 27.4d+4 

c 

 go to (400,400,400,400,400,600,400,600,400,400,400)icm 

c 

 400 continue 

c generic Zry properties 

 ctmelt = 2098.15d0 

 chefus = 22.5d+04 

 go to 50 

 600 continue 

c Zr-1%Nb properties from RRC-KI 

 ctmelt = 2133.d0 

 chefus = 21.0d+4 

 50 continue 

c 

c 5/7/03: ctranb, ctrane, ctranz phase transition temperatures 

c not currently used in FRAPTRAN; RRC-KI does have different 

c transition temperatures for Zr-1%Nb 

 wfox = deloxy + 0.0012d0 

 ctranb = 1094.d0 + wfox*(-1.289d+03 + wfox*7.914d+05) 

 if(wfox .lt. 0.025d0) go to 30 

 ctranb = 1556.4d0 + 3.8281d+04*(wfox - 0.025d0) 

 30 continue 

 ctrane = 392.46d0*((100.d0*deloxy + 0.1242807d0)**2 + 3.1417d0) 

 if(deloxy .lt. 4.7308937d-03) go to 40 

 ctrane = (100.d0*deloxy + 0.12d0)*491.157d0 + 1081.7413d0 

 40 continue 

 ctranz = 1135.15d0 

c 

 return 
A.4 

c 

A.1.2 Fuel Specific Heat Capacity (FCP) and Fuel Enthalpy (FENTHL) 

FRAPCON-3.4: FCP

*deck fcp 

 function fcp (ftemp,facmot,fotmtl,gadoln) 

c 

c the function fcp is used to calculate the specific heat capacity 

c of uo2, puo2, and (u,pu)o2 fuels as a function of temperature, 

c fraction of fuel which is molten, puo2 content, and oxygen-toc metal ratio. 

c fcp = output fuel specific heat capacity (j/(kg*K)) 

c ftemp = input fuel meshpoint temperature (K) 

c facmot = input fuel fraction molten (unitless) 

c facmot = 1.0 - fuel is all molten 

c facmot = 0.0 - fuel is all solid 

c fotmtl = input oxygen to metal ratio (unitless) 

c fotmtl = (atoms oxygen)/(atoms metal) if not known, 

c enter 2.0 

c gadoln = weight fraction of gadolinia in the fuel 

c ufcp = estimated standard error for uo2 (j/(kg*K)) 

c (not currently returned) 

c pufcp = estimated standard error for puo2 (j/(kg*K)) 

c (not currently returned) 

c the equations used in this function are based on data from; 

c (1) j.kerrisk and d.clifton, nuclear technology,16 (1972). 

c (2) o.kruger and h.savage, journal of chemical physics,45 

c (1968). 

c the effect of oxygen-to-metal ratio was estimated from 

c equations published by; 

c (3) c.affortit and j.marcon, revue internationale des hautes 

c temperatures et des refractaires,7 (1970). 

c the specific heat capacity of molten uo2 was taken from 

c (4) l.leibowitz et al, journal of nuclear materials,39 (1970) 

c fcp was originally coded by v.f.baston in march 1974. 

c last modified by g.a.reymann in may 1978. 

c Modification for gadolinia additions were made by DD Lanning 

c in 1996. 

 common / phypro / ftmelt,fhefus,ctmelt,chefus,ctranb, 

 + ctrane,ctranz,fdelta,bu,comp,deloxy 

c the phypro common block is set by call to the subroutine 

c phypro which is part of the materials properties package 

c used at the inel. 

c quantities contained in it used in this subcode are; 

c ftmelt = fuel melting temperature (K) 

c fdelta = liquid-solid coexistence temperature (K) 

c comp = puo2 content (wt.%) 

c when the function fcp is used alone, not as part of the 

c inel code, values for ftmelt, fdelta, and comp must be 

c input. 

 data c1u ,c2u ,c3u ,thu ,edu / 296.7 , 2.43e-02, 8.745e07, 

 + 535.285, 1.577e05 / 

 data c1pu,c2pu,c3pu,thpu,edpu/ 347.4 , 3.95e-04, 3.860e07, 

 + 571.0 , 1.967e05 / 

 data c1gd,c2gd,c3gd,thgd,edgd/ 315.86 , 4.044e-02 , 0.0, 

 + 348.0 , 0.0 / 

 cp(c1,c2,c3,th,ed,t,fotmtl) = c1*(th**2)*exp(th/t)/((t**2)*((exp(t 

 +h/t)-1.)**2))+c2*t+(fotmtl/2.)*(c3*ed/(8.314*(t**2)))*exp(-ed/ 

 +(8.314*t)) 

 t = ftemp 
A.5 

 r = facmot 

 tm = ftmelt 

 fcomp = comp/100.0 

 fcpmol = 503.0 

c fcpmol = specific heat capacity of molten fuel (j/(kg*K)) 

 if (t.gt.(tm+fdelta)) go to 100 

 fcp = cp(c1u,c2u,c3u,thu,edu,t,fotmtl)*(1.0-fcomp)+cp(c1pu,c2pu,c3 

 +pu,thpu,edpu,t,fotmtl)*fcomp 

 fcp = fcp*(1.0-gadoln) + gadoln*cp(c1gd,c2gd,c3gd,thgd,edgd,t, 

 + fotmtl) 

 if (t.lt.(tm-0.1)) go to 110 

 fcp = (1.0-r)*fcp+r*fcpmol 

 go to 110 

100 fcp = fcpmol 

 ufcp = 2.0 

 pufcp = 5.6 

110 continue 

 return 

 end 

FRAPTRAN 1.4: FCP

 

*deck fcp 

 function fcp (ftemp,facmot,fotmtl,frden) 

 implicit real (a-h,o-z) 

c 

c This is a MATPRO-11, Rev. 2 routine modified by PNNL for use in FRAPT 

c 

c the function fcp is used to calculate the specific heat capacity 

c of uo2, puo2, and (u,pu)o2 fuels as a function of temperature, 

c fraction of fuel which is molten, puo2 content, and oxygen-toc metal ratio. 

c 

c fcp = output fuel specific heat capacity (j/(kg*k)) 

c 

c ftemp = input fuel meshpoint temperature (k) 

c facmot = input fuel fraction molten (unitless) 

c facmot = 1.0 - fuel is all molten 

c facmot = 0.0 - fuel is all solid 

c fotmtl = input oxygen to metal ratio (unitless) 

c fotmtl = (atoms oxygen)/(atoms metal) if not known, 

c enter 2.0 

c ufcp = estimated standard error for uo2 (j/(kg*k)) 

c (not currently returned) 

c pufcp = estimated standard error for puo2 (j/(kg*k)) 

c (not currently returned) 

c 

c the equations used in this function are based on data from; 

c (1) j.kerrisk and d.clifton, nuclear technology,16 (1972). 

c (2) o.kruger and h.savage, journal of chemical physics,45 

c (1968). 

c 

c the effect of oxygen-to-metal ratio was estimated from 

c equations published by; 

c (3) c.affortit and j.marcon, revue internationale des hautes 

c temperatures et des refractaires,7 (1970). 

c the specific heat capacity of molten uo2 was taken from 

c (4) l.leibowitz et al, journal of nuclear materials,39 (1970) 

c 

c fcp was originally coded by v.f.baston in march 1974. 

c last modified by g.a.reymann in may 1978. 

c Modified by PNNL, January 1997, to clean up coding and delete lic 
A.6 

c analysis and sensitivity uncertainty coding 

c 

c 

 include 'phypro.h' 

c 

c the phypro common block is set by call to the subroutine 

c phypro which is part of the materials properties package 

c used at the inel. 

c quantities contained in it used in this subcode are; 

c ftmelt = fuel melting temperature (k) 

c fdelta = liquid-solid coexistence temperature (k) 

c compmt = puo2 content (wt.%) 

c when the function fcp is used alone, not as part of the 

c inel code, values for ftmelt, fdelta, and compmt must be 

c input. 

c 

 data c1u ,c2u ,c3u ,thu ,edu / 296.7d0 , 2.43d-02, 8.745d07, 

 & 535.285d0, 1.577d05 / 

 data c1pu,c2pu,c3pu,thpu,edpu/ 347.4d0 , 3.95d-04, 3.860d07, 

 & 571.0d0 , 1.967d05 / 

c 

 cp(c1,c2,c3,th,ed,t,fotmtl) = 

 & c1*(th**2)*exp(th/t)/((t**2)*((exp(th/t)-1.d0)**2)) + c2*t 

 & + (fotmtl/2.d0)*(c3*ed/(8.314d0*(t**2)))*exp(-ed/(8.314d0*t)) 

c 

 t = ftemp 

 r = facmot 

 tm = ftmelt 

 fcomp = compmt/100.0d0 

 fcpmol = 503.0d0 

c fcpmol = specific heat capacity of molten fuel (j/(kg*k)) 

c 

 if (t .gt. (tm+fdelta)) go to 50 

c 

 fcp = cp(c1u ,c2u ,c3u ,thu ,edu ,t,fotmtl)*(1.0d0 - fcomp) + 

 & cp(c1pu,c2pu,c3pu,thpu,edpu,t,fotmtl)*fcomp 

c 

 if (t .lt. (tm-0.1d0)) go to 100 

 fcp = (1.0d0 - r)*fcp + r*fcpmol 

 go to 100 

c 

 50 continue 

 fcp = fcpmol 

 ufcp = 2.0d0 

 pufcp = 5.6d0 

c 

 100 continue 

c 

 return 

c 

 end 

FRAPCON-3.4: FENTHL

*deck fenthl 

c 

 function fenthl (temp,fcomp,fotmtl,ftmelt,facmot,fhefus,gadoln) 

c implicit real * 8 (a-h,o-z) 

c ******************************************************************** 

c this function is called by subroutine energy and computes the 

c enthalpy of fuel at a point relative to zero degrees absolute 

c temperature 
A.7 

c ******************************************************************** 

c input arguments 

c ******************************************************************** 

c facmot - fraction of molten fuel 

c fcomp - puo2 fraction of the fuel 

c fhefus - heat of fussion of the fuel (j/kg) 

c fotmtl - fuel oxygen to metal ratio 

c ftmelt - fuel melting temperature (K) 

c gadoln - weight fraction of gadolinia in the fuel 

c temp - local temperature (K) 

c ******************************************************************** 

c output 

c ******************************************************************** 

c fenthl - local fuel enthalpy relative to zero degrees-K (j/kg) 

c ******************************************************************** 

c fcpmol = specific heat capacity of molten fuel ( j/(kg*K) ) 

 data fcpmol / 503.e0 / 

c the following data statements contain constants from matpro-11 fcp 

 data c1u, c2u, c3u, thu, edu / 296.7 , 2.43e-02, 8.745e07, 

 + 535.285, 1.577e05 / 

 data c1pu,c2pu,c3pu,thpu,edpu / 347.4 , 3.95e-04, 3.860e07, 

 + 571.0 , 1.967e05 / 

 data c1gd,c2gd,c3gd,thgd,edgd / 315.86 , 4.044e-2, 0.0 , 

 + 348.0 , 0.0 / 

c the following equation is the integral of the fuel specific heat 

c with respect to temperature 

 cpdt(c1,th,c2,otm,ed,t,c3) = c1*th*(1./(exp(th/t)-1.))+c2*t*t/2.e0 

 ++c3*otm*exp(-ed/(t*8.314e0))/2.e0 

 tx = temp 

 if (tx.gt.ftmelt) tx = ftmelt 

 fenthl = cpdt(c1u,thu,c2u,fotmtl,edu,tx,c3u)*(1.e0-fcomp) 

 fenthl = fenthl+cpdt(c1pu,thpu,c2pu,fotmtl,edpu,tx,c3pu)*fcomp 

 fenthl = fenthl*(1-gadoln)+ 

 &gadoln*cpdt(c1gd,thgd,c2gd,fotmtl,edgd,tx,c3gd) 

 if (temp.le.ftmelt-2.e0) go to 100 

 fenthl = fenthl+fhefus*facmot 

 if (temp.le.(ftmelt+2.)) go to 100 

 fenthl = fenthl+(temp-ftmelt)*fcpmol 

100 continue 

 return 

 end 

FRAPTRAN 1.4: FENTHL

*deck fenthl 

c 

 function fenthl (temp, fcomp, fotmtl, ftmelt, facmot, fhefus) 

c 

 implicit real (a-h,o-z) 

c 

c 

c ********************************************************************* 

c 

c this function is called by subroutine energy and computes the 

c enthalpy of fuel at a point relative to zero degrees absolute 

c temperature 

c 

c ********************************************************************* 

c input arguments 

c ********************************************************************* 

c 

c facmot - fraction of molten fuel 
A.8 

c fcomp - puo2 fraction of the fuel 

c fhefus - heat of fussion of the fuel (j/kg) 

c fotmtl - fuel oxygen to metal ratio 

c ftmelt - fuel melting temperature (K) 

c temp - local temperature (K) 

c 

c ********************************************************************* 

c output 

c ********************************************************************* 

c 

c fenthl - local fuel enthalpy relative to zero degrees-K (j/kg) 

c 

c ********************************************************************* 

c 

c fcpmol = specific heat capacity of molten fuel ( j/(kg*K) ) 

c 

 data fcpmol / 503.d0 / 

c 

c the following data statements contain constants from matpro-11 fcp 

c 

 data c1u, c2u, c3u, thu, edu / 296.7d0 , 2.43d-02, 8.745d07, 

 & 535.285d0, 1.577d05 / 

 data c1pu,c2pu,c3pu,thpu,edpu / 347.4d0 , 3.95d-04, 3.860d07, 

 & 571.0d0 , 1.967d05 / 

c 

c the following equation is the integral of the fuel specific heat 

c with respect to temperature 

c 

 cpdt(c1,th,c2,otm,ed,t,c3) = 

 & c1*th * ( 1.d0 / (exp(th/t)-1.d0)) 

 & + c2 * t * t / 2.d0 

 & + c3 * otm * exp(-ed /(t * 8.314d0) )/2.d0 

c 

 tx = temp 

 if(tx .gt. ftmelt) tx = ftmelt 

c 

 fenthl = cpdt(c1u ,thu ,c2u ,fotmtl ,edu 

 & ,tx ,c3u ) * ( 1.d0 - fcomp ) 

c 

 if( fcomp .gt. 0.d0 ) fenthl = fenthl + cpdt(c1pu ,thpu ,c2pu 

 & ,fotmtl ,edpu ,tx ,c3pu) * fcomp 

c 

 if( temp .le. ftmelt - 2.d0 ) go to 100 

c 

c 

c 

 fenthl = fenthl + fhefus * facmot 

c 

 if(temp.le.(ftmelt+2.d0))go to 100 

c 

c 

 fenthl = fenthl + (temp - ftmelt) * fcpmol 

 100 continue 

c 

 return 

 end 

A.1.3 Fuel Thermal Conductivity (FTHCON) 

FRAPCON-3.4: FTHCON

*deck fthcon 
A.9 

 subroutine fthcon (ftemp,fraden,fotmtl,con,burnup 

 + ,gadoln,imox) 

c 

c fthcon calculates the fuel thermal conductivity and its 

c derivative with respect to temperature as a function of 

c temperature, density, composition and burnup. 

c 

c UO2 Fuel (IMOX = 0) 

c 

c The equation used in this subroutine is that proposed by 

c staff at NFI, Japan, at the May 1997 ANS Topical Meeting on 

c Light Water Reactor Fuel performance in Portland, OR: (Ohira, 

c K., and N.Itagaki, 1997. "Thermal Conductivity Measurements 

c of High Burnup UO2 Pellet and a Benchmark Calculation of Fuel 

c Center Temperature", proceedings pp. 541-549. Applies to UO2. 

c 

c MOX: 

c 

c Option number 1 (IMOX = 1) 

c 

c The 100% dense solid MOX fuel thermal conductivity formulation is based 

c on a combination of the Duriez stoichiometry-dependent correlation, 

c derived from diffusivity measurements on unirradiated fuel pellets 

c (C.Duriez, et al, J.Nuclear Materials 277, 143-158 2000) and the burnup 

c degradation conatined in a modified version of the NFI fuel thermal 

c conductivity model 

c 

c Option number 2 (IMOX = 2) 

c 

c The MOX fuel thermal conductivity formulation is based 

c on the OECD Halden Reactor Project report "Thermal Performance of 

c of High Burnup Fuel In-pile Temperature Data and Analysis" 

c W.Wiesnack, T. Tverberg, Proceedings of the 2000 International 

c Topical Meeting on LWR Fuel Performance 

c 

c 

c burnup = current local burnup (MWd/MTU) 

c con = output fuel thermal conductivity (W/(m*K)) 

c ftemp = current fuel ring temperature (K) 

c fraden = input fuel density (ratio of actual density to 

c theoretical density) 

c fotmtl = input oxygen to metal ratio of fuel (atoms oxygen/ 

c atoms metal) 

c gadoln = input weight fraction of gadolinia in the fuel 

c 

c the following inputs are by common block 

c comp = input puo2 content of fuel (percent puo2 

c in total fuel weight) 

c bu = input burnup (mw-s/kg-u) 

c emflag(12) = input switch for evaluation model. if this 

c variable is equal to 1.0, the matpro model for 

c fuel thermal conductivity is replaced by the 

c subcode emfton 

c 

 common /uncert/sigftc, sigftex, sigfgr 

 & ,sigswell ,sigcreep ,siggro ,sigcor ,sigh2 

 common / phypro / ftmelt,fhefus,ctmelt,chefus,ctranb, 

 + ctrane,ctranz,fdelta,bu,comp,deloxy 

c 

 include 'lacmdl.h' 

 data on / 1 /, 

 + off / 2 /, 

 + locidx / 12 / 
A.10 

c 

c find constants 

c 

 frpu = comp/100. 

 t = ftemp 

c 

c Burnup in GWd/MTU 

c 

 bug = burnup/1000.0 

c 

 if(imox.eq.0) then 

c 

c NFI formula (Ohira & Itagaki, ANS LWR Fuel perf. Topical mtg. 1997) 

c MODIFIED in January 2002 to raise low-burnup thermal conductivity 

c at low temperature and lower thermal conductivity at very high temp. 

c 

 h = 1/(1.0+396.0*exp(-6380.0/t)) 

 rphonon= 1.0/(0.0452+0.000246*t + 1.0*0.00187*bug+1.1599*gadoln 

 & + (1.0-0.9*exp(-0.04*bug))*0.038*bug**0.28*h) 

 elect = (3.50e9/t**2)*exp(-16361/t) 

 base = rphonon + elect 

c 

c fm is the Lucuta porosity correction factor(applied to 100% TD fuel) 

c 

 fm = fraden/(1.0 + 0.5*(1.0-fraden)) 

c 

c NFI base equation is for 95% TD fuel, so multiply by 1.079 to 

c raise to 100% TD fuel conductivity, then multiply by fm 

c 

 con = base*fm*1.079 

c 

 else if(imox.eq.1) then 

c 

c Using the Duriez/NFI Mod correlation combination 

c 

c base term for MOX 

c where X = deviation from stoichiometry (2-O/M) 

 fm = 1.0789*fraden/(1.0+0.5*(1.0-fraden)) 

c fm is multiplied by 1.0789 to account for 95% TD 

c Porosity correction is Lucuta correction, not Maxwell-Euken 

c as proposed by Duriez et al. 

 x = 2.0-fotmtl 

 ax=2.85*x+0.035 

 cx=(2.86-7.15*x)*1.0e-4 

c 

 h = 1/(1.0+396.0*exp(-6380.0/t)) 

 rphonon = 1.0/(ax + cx*t + 0.00187*bug+1.1599*gadoln 

 &+ (1.0-0.9*exp(-0.04*bug))*0.038*bug**0.28*h) 

 elect = (1.50e9/t**2)*exp(-13520/t) 

 base = rphonon + elect 

 con = base*fm 

c 

 else if(imox.eq.2) then 

c 

c Using the Halden correlation 

c 

 tc=t-273.15 

 tco=min(1650.0,tc) 

 buguo2=bug*0.8815 

 fm = 1.0789*fraden/(1.0+0.5*(1.0-fraden)) 

 base=0.92/(0.1148+0.004*buguo2+1.1599*gadoln+ 

 & 2.475e-4*(1.0-0.00333*buguo2)*tco)+ 

 & 0.0132*exp(0.00188*tc) 
A.11 

 con=base*fm 

c 

c If IMOX.ne.0,1,2 then stop the calculations 

c 

 else 

 stop 'fthcon - IMOX not within bounds' 

 end if 

c 

c find uncertainty 

 if(sigftc.gt.0.0) con = con*(1.0+sigftc*0.088) 

 if(sigftc.lt.0.0) con = con/(1.0-sigftc*0.088) 

 

 if(imox.eq.0) then 

 if(t.lt.ftmelt) then 

 ucon = 0.2*(1.0+abs(2.0-fotmtl)*10.) 

 else 

 ucon = con/2.0 

 endif 

 else 

 if(t.le.1800.0) then 

 ucon = 0.07*con 

 else 

 frac=(t-1800.0)/(3100.0-1800.0)*(0.20-0.07)+0.07 

 ucon=frac*con 

 endif 

 endif 

 if (emflag(locidx).eq.on) call emfton (ftemp,fraden,ftmelt,con) 

 return 

 end 

FRAPTRAN 1.4: FTHCON

*deck fthcon 

c 

 subroutine fthcon (ftemp, fraden, fotmtl, con, dkdt, burnup, 

 & gadoln, compmt) 

c 

 implicit real (a-h, o-z) 

c 

c FTHCON calculates the UO2 fuel thermal conductivity 

c as a function of temperature, density, composition and burnup. 

c 

c The equation used in this subroutine is that proposed by 

c staff at NFI, Japan, at the May 1997 ANS Topical Meeting on 

c Light Water Reactor Fuel Performance in Portland, Oregon. 

c reference: Ohira and Itakaki, 1997. "Thermal Conductivity 

c Measurements of High Burnup UO2 Pellet and a Benchmark 

c Calculation of Fuel Center Temperature," proceedings pg. 541-549. 

c 

c The NFI formula was modified in January 2002 by PNNL to raise 

c low-burnup thermal conductivity at low temperature and lower thermal 

c conductivity at very high temperature. 

c The PNNL reference is: Lanning and Beyer. 2002. "Revised UO2 Thermal 

c Conductivity for FRAPCON-3 NRC Fuel Performance Codes," in Proceedings of 

c ANS Annual Meeting, June 9-13, 2002, Hollywood, Florida. 

c 

c This is the same formulation used in FRAPCON-3.2 (February 2003), 

c minus the formulations for MOX fuel. 

c 

c burnup = current local burnup (MWd/MTU) 

c con = output fue thermal conductivity (W/(m*K)) 
A.12 

c ftemp = current fuel ring temperature (K) 

c fraden = input fuel density (ratio of actual to theoretical) 

c fotmtl = input oxygen to metal ratio of fuel 

c (atoms oxygen/atoms metal) 

c gadoln = input weight fraction of gadolinia in the fuel 

c compmt = input weight % of Pu in Pu and U 

c 

c convert burnup to GWd/MTU 

 bug = burnup/1.d3 

 if(bug .lt. 1.d-10) bug = 0.001 

c use modified NFI model for UO2 and Duriez/NFI Mod model for MOX 

 if(compmt.gt.0.0) go to 100 

c UO2 Model 

c 

c NFI formula (Ohira & Itagaki, ANS LWR Fuel perf. Topical mtg. 1997) 

c MODIFIED in January 2002 to raise low-burnup thermal conductivity 

c at low temperature and lower thermal conductivity at very high temp. 

c 

c calculate terms 

c temperature dependence of annealing of irradiation effects 

 h = 1.0/(1.0 + 396.0*exp(-6380.0/ftemp)) 

c phonon term 

 rphonon = 1.0/(0.04520 + 0.000246*ftemp + 0.00187*bug 

 & + (1.0 - 0.9*exp(-0.04*bug))*(0.038*bug**0.28)*h) 

c electronic term 

 elect = (3.50d9/ftemp**2)*exp(-16361./ftemp) 

c base is conductivity before apply porosity correction 

 base = rphonon + elect 

c apply lucuta porosity correction factor (applied to 100% TD fuel) 

 fm = fraden/(1.0 + 0.5*(1.0-fraden)) 

c NFI base equation is for 95% TD fuel, so multiply by 1.079 to 

c raise to 100% TD fuel conductivity, then multiply by fm 

 con = base*fm*1.079 

 go to 200 

c 

c MOX Model 

c 

100 continue 

c Using the Duriez/NFI Mod correlation combination 

c 

c base term for MOX 

c where X = deviation from stoichiometry (2-O/M) 

 fm = 1.0789*fraden/(1.0+0.5*(1.0-fraden)) 

c fm is multiplied by 1.0789 to account for 95% TD 

c Porosity correction is Lucuta correction, not Maxwell-Euken 

c as proposed by Duriez et al. 

 x = 2.0-fotmtl 

 ax=2.85*x+0.035 

 cx=(2.86-7.15*x)*1.0e-4 

c 

 h = 1/(1.0+396.0*exp(-6380.0/ftemp)) 

 rphonon = 1.0/(ax + cx*ftemp + 0.00187*bug+1.1599*gadoln 

 &+ (1.0-0.9*exp(-0.04*bug))*0.038*bug**0.28*h) 

 elect = (1.50e9/ftemp**2)*exp(-13520/ftemp) 

 base = rphonon + elect 

 con = base*fm 

c 

200 continue 

 return 

 end 

A.1.4 Fuel Emissivity (FEMISS) 
A.13 

FRAPCON-3.4: FEMISS

*deck femiss 

 function femiss (ftemp) 

c 

c femiss calculates fuel emissivity as a function of temperature. 

c ftemp = input fuel temperature (k) 

c femiss = output fuel emissivity (unitless) 

c data used to develop the model are 

c held and wilder, journ. amer. ceram. soc. 

c vol. 52, (1969) 

c cabannes, et at, c. r. acad. sci., paris, ser. b 

c (1967) 

c femiss was coded by r. e. mason in october 1978. 

 femiss = 0.78557+1.5263e-05*ftemp 

c the following calculations provide the upper and lower bounds. 

c the upper and lower bounds are not an output unless the user 

c desires to modify the subcode appropriately. 

c femisu is the upper bound 

 femisu = femiss*(1.+0.06796) 

c femisl is the lower bound 

 femisl = femiss*(1.-0.06796) 

100 continue 

 return 

 end 

FRAPTRAN 1.4: FEMISS

*deck femiss 

 function femiss (ftemp) 

 implicit real (a-h,o-z) 

c 

c This is a MATPRO-11, Rev. 2 routine modified by PNNL for use in FRAPT 

c 

c femiss calculates fuel emissivity as a function of temperature. 

c 

c ftemp = input fuel temperature (k) 

c femiss = output fuel emissivity (unitless) 

c 

c data used to develop the model are 

c held and wilder, journ. amer. ceram. soc. 

c vol. 52, (1969) 

c cabannes, et at, c. r. acad. sci., paris, ser. b 

c (1967) 

c 

c femiss was coded by r. e. mason in october 1978. 

c Modified by PNNL, January 1997, to clean up coding and delete lic 

c analysis and sensitivity uncertainty coding 

c 

 femiss = 0.78557d0 + 1.5263d-05 * ftemp 

c 

c the following calculations provide the upper and lower bounds. 

c the upper and lower bounds are not an output unless the user 

c desires to modify the subcode appropriately. 

c femisu is the upper bound 

 femisu = femiss * (1.d0 + 0.06796d0) 

c femisl is the lower bound 

 femisl = femiss * (1.d0 - 0.06796d0) 

c 

 return 

c 
A.14 

 end 

A.1.5 Fuel Thermal Expansion (FTHEXP) 

FRAPCON-3.4: FTHEXP

*deck fthexp 

 function fthexp (ftemp,facmot) 

c 

c the subprogram fthexp is used to calculate the strain of the 

c fuel caused by thermal expansion. 

c thermal expansion strain is calculated for uo2,pu02, or (u,pu)02 

c as a function of temperature, puo2 content, and fraction of 

c the fuel which is molten. 

c the correlations were derived with the assumption that the 

c thermal expansion strain is zero at 300k. 

c fthexp = output fuel strain due to thermal expansion (unitless) 

c ufthex = estimated standard error of fthexp (per cent) 

c (not currently returned). 

c ftemp = input fuel temperature (k) 

c facmot = input fuel fraction which is molten (unitless) 

c facmot = 1.0 - fuel is all molten 

c facmot = 0.0 - fuel is all solid 

c the correlations used in this subroutine are based on data from; 

c (1) p.j.baldock et al, journal of nuclear materials, 18 (1966) 

c (2) n.h.brett and l.e.russel, proceedings of the 2nd 

c international conference on plutonium metallurgy, grenoble, 

c france (1960) pp 397-410. 

c (3) m.d.burdock and h.s.parker, journal of the american 

c ceramic society,39 (1956) pp 181-187. 

c (4) j.a.christensen, journal of the american ceramic society, 

c 46 (1963) pp 607-608. 

c (5) j.b.conway et al, transactions of the american nuclear 

c society,6 (1963). 

c (6) f.gronvold, journal of inorganic and nuclear chemistry,1 

c (1955) pp 357-370. 

c (7) m.hoch and a.c.momin, high temperatures-high pressures,1 

c (1969) pp 401-407. 

c (8) c.p.kempter and r.o.elliott, the journal of chemical 

c physics,30 (1958) pp 1524-1526. 

c (9) w.a.lambertson and j.h.handwerk, anl-5053 (1956). 

c (10)m.tokar et al, nuclear technology, 17 (1973) pp 147-152. 

c fthexp was originally coded by v.f.baston in march 1974 

c modified by c.s.olsen in feb. 1975 

c last modified by g.a.reymann in july 1978 

 common /uncert/sigftc, sigftex, sigfgr 

 & ,sigswell ,sigcreep ,siggro ,sigcor ,sigh2 

 common / phypro / ftmelt,fhefus,ctmelt,chefus,ctranb, 

 + ctrane,ctranz,fdelta,bu,comp,deloxy 

c the phypro common block is set by call to the subroutine 

c phypro which is part of the materials properties package 

c used at the inel. 

c quantities contained in it used in this subcode are; 

c comp = puo2 content (wt. :) 

c ftmelt = fuel melting temperature (k) 

c fdelta = liquid-solid coexistence temperature range (k) 

c when the function fthexp is used alone, not as part of 

c the inel code, values for comp, ftmelt, and fdelta must 

c be input. 

c 

c c1u, c2u, c3u, edu parameters updated by W.G. Luscher 10-17-07 

c 
A.15 

 data c1u ,c2u ,c3u ,edu / 9.8e-06, 2.61e-03, 3.16e-01, 1.32e-19 / 

 data c1pu,c2pu,c3pu,edpu/ 9.0e-06, 2.7e-03, 7.0e-02, 7.0e-20 / 

 ftx(c1,c2,c3,ed,bk,t) = c1*t-c2+c3*exp(-ed/(bk*t)) 

 t = ftemp 

 bk = 1.38e-23 

c bk = boltzmann's constant (j/k) 

 fcomp = comp/100.0 

 if (t.gt.ftmelt) go to 100 

 fthexp = ftx(c1u,c2u,c3u,edu,bk,t)*(1.0-fcomp)+ftx(c1pu,c2pu,c3pu, 

 +edpu,bk,t)*fcomp 

 go to 120 

100 fthexm = ftx(c1u,c2u,c3u,edu,bk,ftmelt)*(1.0-fcomp)+ftx(c1pu,c2pu, 

 +c3pu,edpu,bk,ftmelt)*(fcomp) 

 if (t.ge.(ftmelt+fdelta)) go to 110 

 fthexp = fthexm+.043*facmot 

 go to 120 

110 fthexp = fthexm+.043+3.6e-05*(t-(ftmelt+fdelta)) 

120 continue 

c add on uncertainty 

 if(sigftex.gt.0.0) fthexp = fthexp*(1.0+sigftex*0.103) 

 if(sigftex.lt.0.0) fthexp = fthexp/(1.0-sigftex*0.103) 

 ufthex = 10.0 

 return 

 end 

FRAPTRAN 1.4: FTHEXP

*deck fthexp 

c 

 function fthexp (ftemp,facmot) 

c 

 implicit real (a-h,o-z) 

c 

c This is a MATPRO-11, Rev. 2 routine modified by PNNL for use in FRAPT 

c 

c the subprogram fthexp is used to calculate the strain of the 

c fuel caused by thermal expansion. 

c thermal expansion strain is calculated for uo2,pu02, or (u,pu)02 

c as a function of temperature, puo2 content, and fraction of 

c the fuel which is molten. 

c the correlations were derived with the assumption that the 

c thermal expansion strain is zero at 300k. 

c 

c fthexp = output fuel strain due to thermal expansion (unitless) 

c ftemp = input fuel temperature (k) 

c facmot = input fuel fraction which is molten (unitless) 

c facmot = 1.0 - fuel is all molten 

c facmot = 0.0 - fuel is all solid 

c 

c the correlations used in this subroutine are based on data from; 

c (1) p.j.baldock et al, journal of nuclear materials, 18 (1966) 

c (2) n.h.brett and l.e.russel, proceedings of the 2nd 

c international conference on plutonium metallurgy, grenoble, 

c france (1960) pp 397-410. 

c (3) m.d.burdock and h.s.parker, journal of the american 

c ceramic society,39 (1956) pp 181-187. 

c (4) j.a.christensen, journal of the american ceramic society, 

c 46 (1963) pp 607-608. 

c (5) j.b.conway et al, transactions of the american nuclear 

c society,6 (1963). 

c (6) f.gronvold, journal of inorganic and nuclear chemistry,1 

c (1955) pp 357-370. 

c (7) m.hoch and a.c.momin, high temperatures-high pressures,1 
A.16 

c (1969) pp 401-407. 

c (8) c.p.kempter and r.o.elliott, the journal of chemical 

c physics,30 (1958) pp 1524-1526. 

c (9) w.a.lambertson and j.h.handwerk, anl-5053 (1956). 

c (10)m.tokar et al, nuclear technology, 17 (1973) pp 147-152. 

c 

c fthexp was originally coded by v.f.baston in march 1974 

c modified by c.s.olsen in feb. 1975 

c last modified by g.a.reymann in july 1978 

c Modified by PNNL, January 1997, to clean up coding and delete lic 

c analysis and sensitivity uncertainty coding 

c 

 include 'phypro.h' 

c 

c the phypro common block is set by call to the subroutine 

c phypro which is part of the materials properties package 

c used at the inel. 

c quantities contained in it used in this subcode are; 

c compmt = puo2 content (wt. :) 

c ftmelt = fuel melting temperature (k) 

c fdelta = liquid-solid coexistence temperature range (k) 

c when the function fthexp is used alone, not as part of 

c the inel code, values for compmt, ftmelt, and fdelta must 

c be input. 

c 

c c1u, c2u, c3u, and edu parameters updated by W.G. Luscher 10-17-07 

c 

 data c1u ,c2u ,c3u ,edu / 9.8d-06, 2.61d-03, 3.16d-01, 1.32d-19 / 

 data c1pu,c2pu,c3pu,edpu/ 9.0d-06, 2.7d-03, 7.0d-02, 7.0d-20 / 

c bk is Boltzmann's constant (J/k) 

 data bk / 1.38d-23 / 

c 

 ftx (c1,c2,c3,ed,bk,t) = c1 * t - c2 + c3 * exp(-ed/(bk * t)) 

c 

 t = ftemp 

 fcomp = compmt/100.0d0 

c 

 if(t .gt. ftmelt) go to 25 

c 

 fthexp = ftx(c1u ,c2u ,c3u ,edu ,bk,t) * (1.0d0 - fcomp) 

 & + ftx(c1pu,c2pu,c3pu,edpu,bk,t) * fcomp 

 go to 100 

c 

 25 continue 

 fthexm = ftx(c1u ,c2u ,c3u ,edu ,bk,ftmelt) * (1.0d0 - fcomp) 

 & + ftx(c1pu,c2pu,c3pu,edpu,bk,ftmelt) * fcomp 

c 

 if(t .ge. (ftmelt+fdelta)) go to 75 

c 

 fthexp = fthexm + 0.043d0 * facmot 

 go to 100 

c 

 75 continue 

 fthexp = fthexm + 0.043d0 + 3.6d-05 * (t- (ftmelt + fdelta)) 

c 

 100 continue 

c 

 return 

 end 

A.1.6 Fuel Densification (FUDENS) 
A.17 

FRAPCON-3.4: FUDENS

*deck fudens 

 function fudens (ftemp,bu,fdens,rsntr,tsint,comp,prvden) 

c 

c fudens calculates irradiation-induced densification. 

c fudens = output fuel dimensional change (%) 

c ftemp = input fuel temperature (k) 

c bu = input burnup (mw-s/kg-u) 

c fdens = input fuel density (kg/m**3) 

c rsntr = input maximum density change determined by a resintering 

c test of 1973 k for 24 hours (kg/m**3) 

c tsint = input fuel sintering temperature (k) 

c comp = input plutonia content (weight percent) 

c prvden = input total densification from previous time step (%) 

c fudens was developed and programmed by c. s. olsen (january 1975). 

c updated and corrected by b. w. burnham (october 1975). 

c fudens was modified by r. e. mason (november 1978). 

 dimension c(2),b(5) 

 data c / 10.96, 11.46/ 

 data b / 3.0, 1.00, 3.0, 2.00, 35.00 / 

 dlen2(alen,bu,abu) = -b(1)+alen+b(2)*exp(-b(3)*(bu+abu))+b(4)*exp( 

 +-b(5)*(bu+abu)) 

 dlen3(bu) = -b(2)*b(3)*exp(-b(3)*bu)-b(4)*b(5)*exp(-b(5)*bu) 

 fbu = bu*1.157e-05 

 ts = tsint-2.7315e02 

c if rsntr or tsint is not defined by user, the default value is 

c tsint = 1873 k. 

 if (tsint.le.0.0) ts = 1600.0 

 roth = c(1)*c(2)/(0.01*comp*c(1)+(1.0-0.01*comp)*c(2)) 

 de = fdens/(roth*10.0) 

 if (rsntr .gt. 0.0) dlen1 = 100.0 * rsntr/(3.0 * fdens) 

 if ((ftemp.ge.1000.).and.(rsntr.le.0.)) dlen1 = 66.6*(100.0-de)/(t 

 +s-1180.0) 

 if ((ftemp.lt.1000.).and.(rsntr.le.0.)) dlen1 = 22.2*(100.0-de)/(t 

 +s-1180.0) 

 x3 = 0.0 

 x4 = 1.0 

 al1 = dlen1 

 al3 = 3.0-al1 

 al4 = 0.0 

 if (al3.le.4.27e-03) go to 130 

 do 100 i=1,6 

 y2 = dlen2(al3,x4,al4) 

 y1 = dlen2(al3,x3,al4) 

 if (y1*y2.le.0.) go to 110 

 x3 = x4 

 x4 = x4+1.0 

 if (i.eq.6) go to 160 

100 continue 

110 continue 

 x1 = x3 

 do 120 j=1,50 

 x = x1-dlen2(al3,x1,al4)/dlen3(x1) 

 err = abs((x-x1)*100.0/x) 

 if (err.le.2.0e-04) go to 140 

 x1 = x 

120 continue 

130 al3 = 2.996 

 al2 = 5.384 

 go to 150 

140 al2 = x 

150 continue 
A.18 

 fuden = dlen2(al3,fbu,al2) 

 if (bu.lt.1728) fuden = 0.0 

 go to 170 

160 print 180 

 fuden = 0.0 

170 continue 

 if (abs(fuden).le.abs(prvden)) fudens = 0.0 

 if (abs(fuden).gt.abs(prvden)) fudens = fuden-prvden 

 return 

180 format (1x,/' no roots found between 0 and 6000 mws/mt uo2') 

 end 

A.1.7 Fuel Swelling (FSWELL) 

FRAPCON-3.4: FSWELL

*deck fswell 

 subroutine fswell (fdens,comp,bu,bul,ftemp,soldsw) 

c 

c comp = input plutonia content (weight percent) 

c fdens = input initial density of the fuel (kg/m**3) 

c ftemp = input temperature of the fuel ring (k) 

c bu = input burnup to end of time step (mw-s/kg-u) 

c bul = input burnup to end of last time step (mw-s/kg-u) 

c soldsw = output fuel swelling due to solid fission products 

c (fraction) 

c fswell was developed and programmed by r. e. mason - june 1978. 

c fswell was revised July 1995 by G. A. Berna per D. D. Lanning 

c recommendations in ltr to L. Siefken, May 4, 1995 

c fswell was revised October 2008 by K.J. Geelhood 

c - program - 

 common /uncert/sigftc, sigftex, sigfgr 

 & ,sigswell ,sigcreep ,siggro ,sigcor ,sigh2 

 bus = fdens*2.974e+10*(bu-bul) 

 soldsw = bus*(2.315e-23+sigswell*2.315e-24) 

 if(bu.ge.6912000.0) soldsw = bus*(3.211e-23+sigswell*3.211e-24) 

 return 

 end 

A.2 Cladding Material Properties 

A.2.1 Cladding Specific Heat (CCP) 

FRAPCON-3.4: CCP

*deck ccp 

 function ccp (ctemp) 

c 

c ccp calculates the specific heat at constant pressure 

c for zircaloys 

c ccp = output cladding specific heat at constant 

c pressure (j/kg-k) 

c ctemp = input cladding temperature (k) 

c this code is based on data from 

c (1) c. r. brooks and e. e. stansbury, "the specific heat 0f 

c of zircaloy-2 from 50 to 700 c", journal of nuclear 

c materials 18 (1966) p 233 

c (2) e. a. eldridge and h. w. deem, specific heats and heats 

c of transformation of zircaloy-2 and low nickel zircaloy-2 
A.19 

c bmi-1803 (may 31, 1967) 

c conversion from j/(kg*k) to cal/(gm*c) is 

c 2.39006d-4 (cal/(gm*c))/(j/(kg*k)) 

c estimated standard error of the prediction of 

c ccp for the specific heat of zircaloy cladding 

c samples -- (sum of squared residuals/(number of 

c residuals-degrees of freedom))**0.5 -- is 

c (1) for temperature less than 1090 k, 10j/kg-k) 

c (2) for temperature 1090 k to 1300 k, 25j/kg-k) 

c (3) for temperature above 1300 k , 100j/kg-k) 

c ccp coded by r. l. miller october 1974 

c modified by d. l. hagrman may 1976 

 include 'lacmdl.h' 

 dimension cpdata(26) 

 data cpdata/ 281., 300., 302., 400., 331.,640., 

 + 375.,1090., 502.,1093., 590.,1113., 615.,1133., 

 + 719.,1153., 816.,1173., 770.,1193., 619.,1213., 

 + 469.,1233., 356.,1248. / 

 data npcp, iu / 13,1/ 

 data on / 1 /, 

 + off / 2 /, 

 + locidx / 4 / 

 if (emflag(locidx).eq.on) go to 110 

 if (ctemp.ge.1248.0) go to 100 

 ccp = polate(cpdata,ctemp,npcp,iu) 

 go to 120 

100 ccp = 356. 

 go to 120 

110 ccp = emccp(ctemp) 

120 continue 

 return 

 end 

FRAPTRAN 1.4: CCP

*deck ccp 

 function ccp (ctemp, icm) 

 implicit real (a-h,o-z) 

c 

c This is a MATPRO-11, Rev. 2 routine modified by PNNL for use in FRAP 

c 

c ccp calculates the specific heat at constant pressure 

c for zircaloys 

c 

c ccp = output cladding specific heat at constant 

c pressure (j/kg-k) 

c ctemp = input cladding temperature (k) 

c 

c this code is based on data from 

c (1) c. r. brooks and e. e. stansbury, "the specific heat 0f 

c of zircaloy-2 from 50 to 700 c", journal of nuclear 

c materials 18 (1966) p 233 

c (2) e. a. eldridge and h. w. deem, specific heats and heats 

c of transformation of zircaloy-2 and low nickel zircaloy-2 

c bmi-1803 (may 31, 1967) 

c 

c conversion from j/(kg*k) to cal/(gm*c) is 

c 2.39006d-4 (cal/(gm*c))/(j/(kg*k)) 

c 

c estimated standard error of the prediction of 

c ccp for the specific heat of zircaloy cladding 

c samples -- (sum of squared residuals/(number of 
A.20 

c residuals-degrees of freedom))**0.5 -- is 

c (1) for temperature less than 1090 k, 10j/kg-k) 

c (2) for temperature 1090 k to 1300 k, 25j/kg-k) 

c (3) for temperature above 1300 k , 100j/kg-k) 

c 

c ccp coded by r. l. miller october 1974 

c modified by d. l. hagrman may 1976 

c Modified by PNNL, January 1997, to clean up coding, including 

c removed licensing assistance coding and sensitivity uncertainty 

c coding 

c 

c Zry properties 

 dimension cpdata(28) 

 data cpdata/279.d0,290.d0, 

 & 281.d0, 300.d0, 302.d0, 400.d0, 331.d0, 640.d0, 

 & 375.d0,1090.d0, 502.d0,1093.d0, 590.d0,1113.d0, 615.d0,1133.d0, 

 & 719.d0,1153.d0, 816.d0,1173.d0, 770.d0,1193.d0, 619.d0,1213.d0, 

 & 469.d0,1233.d0, 356.d0,1248.d0 / 

 data npcp, iu / 14 , 1 / 

c 

c Zr-1%Nb properties 

c low heat rate (<1000K/s) 

 dimension cpwwer(24) 

 data cpwwer/345.0d0,280.d0, 

 & 360.d0, 473.d0, 370.d0, 573.d0, 380.d0, 673.d0, 

 & 383.d0, 773.d0, 385.d0, 873.d0, 448.d0, 883.d0, 680.d0, 973.d0, 

 & 816.d0,1025.d0, 770.d0,1073.d0, 400.d0,1153.d0, 392.d0,1173.d0/ 

 data npww / 12 / 

c high heat rate (>1000K/s) 

 dimension cpwwr2(26) 

 data cpwwr2/412.6d0,1100.d0, 

 & 420.d0,1110.d0, 480.d0,1120.d0, 600.d0,1134.d0, 

 & 1000.d0,1142.d0, 1400.d0,1150.d0,1600.d0,1155.d0,1400.d0,1161.d0, 

 & 1000.d0,1168.d0, 600.d0,1177.d0, 400.d0,1180.d0, 360.d0,1190.d0, 

 & 348.d0,1200.d0/ 

 data npw2 / 13 / 

c 

 go to (400,400,400,400,400,600,400,600,400,400,400)icm 

c 

 400 continue 

c generic Zry 

 ccp = 279.d0 

 if (ctemp.ge.290.d0 .and. ctemp.le.1248.d0) 

 & ccp = polate (cpdata, ctemp, npcp, iu) 

 if (ctemp .ge. 1248.d0) ccp = 356.d0 

 return 

c 

 600 continue 

c RRC-KI Zr-1%Nb property 

c 

c This block calculates the specific heat at constant pressure for 

c H1-alloy cladding. If iwwer = 1 then data is derived from 

c Volkov B.Yu. et. al."Material property library for H1-alloy 

c cladding", Preprint IAE-4941/11, 1989. This property are induced 

c by low heating rate. If iwwer = 2 then data is derived from 

c Ljusternik V.E. et. al."Experimental research of zirconium reactor 

c materials thermal properties:110-alloy". M.,J. High Temperature 

c Thermal Physics." v.31 n.4, 1993. This property are induced by 

c high-rate heating>=1000 K/s. 

c 

c 5/5/03: default to iwwer=1 until can set this flag up 

 iwwer = 1 

c low heat rate 
A.21 

 if(iwwer.eq.2) go to 601 

 if(ctemp .lt. 1173.d0) ccp = polate(cpwwer,ctemp,npww,iu) 

 if(ctemp.ge.1173.d0) ccp = 392.d0 

 return 

c 

 601 continue 

c high heat rate 

 ccp = 237.5d0+15.91d-2*ctemp 

 if(ctemp.ge.1100.d0 .and. ctemp.le.1200.d0) 

 & ccp = polate(cpwwr2,ctemp,npw2,iu) 

 if(ctemp.gt.1200.d0) ccp = 199.7d0+12.364d-2*ctemp 

 return 

c 

 end 

A.2.2 Cladding Thermal Conductivity (CTHCON) 

FRAPCON-3.4: (CTHCON)

*deck cthcon 

 subroutine cthcon (ctemp,time,flux,coldw,ccon) 

c 

c cthcon calculates cladding thermal conductivity as a function 

c of temperature, time, flux, and cold work 

c ccon = output thermal conductivity of zircaloy-4 (w/m-k) 

c ctemp = input cladding meshpoint temperature (k) 

c time = input time at temperature and flux (s) 

c flux = input fast neutron flux ((neutrons/m**2)/s) 

c coldw = input cold work (unitless ratio of areas) 

c the equation used in this subcode is based on data from 

c w.k. anderson, c.j. beck, a.r. kephart and j.s. theilacker 

c astm-stp-314, 1962, pp 62 - 93 

c j.n. chirigos, c. kass, w.s. kirk and g.j. salvaggio 

c fuel element fabrication, academic press, 1961, pp 19 - 55 

c a.d. feith, gemp-669 (oct 1966) 

c c.f. lucks and h.w. deem, bmi-1273 (1958) pp 7-9 

c a.e. powers, kapl-2146 (1961) 

c d.b. scott, wcap-3269-41 (1965) pp 5-9 

c this version of cthcon does not use time, flux or coldwork 

c as parameters in calculation of zircaloy thermal conductivity 

c one standard deviation of this function = 1.01 w/m-k 

c this version of cthcon was developed by r.l. miller, dec 1975 

c last modified by e.r. carlson in june 1978. 

 include 'lacmdl.h' 

 data on / 1 /, 

 + off / 2 /, 

 + locidx / 7 / 

 if (emflag(locidx).eq.on) go to 100 

 ccon = 7.511+ctemp*(2.088e-2+ctemp*(-1.450e-5+ctemp*7.668e-09)) 

 go to 110 

100 call emcton (ctemp,time,flux,coldw,ccon) 

110 continue 

 return 

 end 

FRAPTRAN 1.4: (CTHCON)

*deck cthcon 

 subroutine cthcon (ctemp, time, flux, coldw, ccon, cdkdt, icm) 

 implicit real (a-h,o-z) 

c 
A.22 

c This is a MATPRO-11, Rev. 2 routine modified by PNNL for use in FRAPT 

c 

c cthcon calculates cladding thermal conductivity as a function 

c of temperature, time, flux, and cold work 

c 

c ccon = output thermal conductivity of zircaloy-4 (w/m-k) 

c cdkdt = output derivative of thermal conductivity with 

c respect to temperature 

c 

c ctemp = input cladding meshpoint temperature (k) 

c time = input time at temperature and flux (s) 

c flux = input fast neutron flux ((neutrons/m**2)/s) 

c coldw = input cold work (unitless ratio of areas) 

c 

c go to appropriate model for cladding type 

 go to (400,400,400,400,400,600,400,600,400,400,400) icm 

c 

 400 continue 

c Generic Zircaloy property from MATPRO 

c 

c the equation used in this subcode is based on data from 

c w.k. anderson, c.j. beck, a.r. kephart and j.s. theilacker 

c astm-stp-314, 1962, pp 62 - 93 

c j.n. chirigos, c. kass, w.s. kirk and g.j. salvaggio 

c fuel element fabrication, academic press, 1961, pp 19 - 55 

c a.d. feith, gemp-669 (oct 1966) 

c c.f. lucks and h.w. deem, bmi-1273 (1958) pp 7-9 

c a.e. powers, kapl-2146 (1961) 

c d.b. scott, wcap-3269-41 (1965) pp 5-9 

c 

c this version of cthcon does not use time, flux or coldwork 

c as parameters in calculation of zircaloy thermal conductivity 

c 

c one standard deviation of this function = 1.01 w/m-k 

c 

c this version of cthcon was developed by r.l. miller, dec 1975 

c last modified by e.r. carlson in june 1978. 

c Modified by PNNL, January 1997, to clean up coding and delete 

c licening analysis and sensitivity uncertainty coding 

c 

 ccon = 7.511d0 + ctemp*(2.088d-2 + ctemp * (-1.450d-5 + 

 & ctemp*7.668d-09)) 

 cdkdt = 2.088d-2 + ctemp*(-2.9d-5 + ctemp*2.3d-8) 

 return 

c 

 600 continue 

c E-110 property as coded by RRC-KI and provided to PNNL 

c 

c This code block calculates thermal conductivity (w/(m*k)) and 

c derivative of thermal conductivity with respect to temperature 

c (w/(m*k*k)) of H1-alloy. Data is derived from: 

c Volkov B.Yu. et. al."Material property library for H1-alloy 

c cladding", Preprint IAE-4941/11, 1989. 

c 

 if(ctemp.le.2133.0d0) then 

 ccon = dexp(0.000461843d0*ctemp)*15.0636d0 

 cdkdt = 0.006957018215d0*dexp(0.000461843d0*ctemp) 

 else 

 ccon = 36.d0 

 cdkdt = 0.d0 

 endif 

 return 

c 
A.23 

 end 

A.2.3 Cladding Oxide Thermal Conductivity (ZOTCON) 

FRAPCON-3.4: ZOTCON

*deck zotcon 

 function zotcon (ctemp) 

c 

c zotcon calculates zirconium dioxide thermal conductivity as a 

c function of temperature. 

c zotcon = output thermal conductivity of zro2 (w/m-k) 

c ctemp = input cladding meshpoint temperature (k) 

c the data used to generate this correlation were taken from: 

c w.d. kingery, j. francl, r.l. coble and t. vasilos 

c j. amer. ceram. soc., 37 no.i (1954) pp 107-110 

c zotcon was coded by r.l.miller in march 1974. 

c modified by r.l. miller may 1974 

c use caution above 1700k 

 t = ctemp 

 zotcon = 1.9599-t*(2.41e-4-t*(6.43e-7-t*1.946e-10)) 

 return 

 end 

FRAPTRAN 1.4: ZOTCON

*deck zotcon 

 function zotcon (ctemp) 

 implicit real (a-h,o-z) 

c 

c zotcon calculates zirconium dioxide thermal conductivity as a 

c function of temperature. 

c 

c zotcon = output thermal conductivity of zro2 (W/m-K) 

c ctemp = input cladding meshpoint temperature (K) 

c 

c the data used to generate this correlation were taken from: 

c w.d. kingery, j. francl, r.l. coble and t. vasilos 

c j. amer. ceram. soc., 37 no.i (1954) pp 107-110 

c 

c zotcon was coded by r.l.miller in march 1974. 

c modified by r.l. miller may 1974 

c 

c use caution above 1700K 

c 

 t = ctemp 

 zotcon = 1.9599d0 - t*(2.41d-4 - t*(6.43d-7 - t*1.946d-10)) 

c 

 return 

 end 

A.2.4 Cladding Oxide Emissivity (ZOEMIS) 

FRAPCON-3.4: ZOEMIS

*deck zoemis 

 subroutine zoemis (ctmax,zroxid,emissv) 

c 

c zoemis calculates the emissivity of the cladding surface 
A.24 

c as a function of maximum cladding temperature and 

c oxide thickness. expected standard error of the 

c calculated emissivity compared to in-reactor 

c data is also returned. 

c emissv = output cladding surface emissivity (unitless) 

c puemis = output positive standard error expected in emissv 

c when compared to in-reactor data (not currently returned) 

c uuemis = output negative standard error expected in emissv 

c when compared to in-reactor data (not currently returned) 

c ctmax = input maximum cladding temperature (k) 

c zroxid = input oxide layer thickness (m) 

c the equations used in this subroutine are based on data from 

c (1) aec fuels and materials development program progress 

c report no. 76, usaec report gemp - 1008 (1968). section by 

c e. f. juenke and s. j. sjodahl, p 239 

c (2) t. b. burgoyne and a. garlick, paper presented at 

c specialists meeting on the behaviour of water reactor 

c fuel elements under accident conditions, spatind norway 

c (september 1976) 

c (3) e. v. murphy and f. havelock, emissivity of zirconium 

c alloys in air in the temperature range 100 - 400 c 

c j. nuc. mat., 60 (1976) pp 167-176 

c zoemis coded by r. l. miller sept 1974 

c modified by d. l. hagrman october 1976 

c model for temperatures below 1500k follows 

 if (zroxid.ge.3.88e-06) go to 100 

 emissv = 3.25e-01+1.246e05*zroxid 

 go to 110 

100 emissv = 8.08642e-01-5.00e01*zroxid 

110 puemis = 0.1 

 uuemis = 0.1 

 if (ctmax.le.1500.) go to 120 

c modification for maximum temperatures above 1500 k follows 

 emissv = emissv*exp((1.50e03-ctmax)/3.00e02) 

 if (emissv.lt.0.325) emissv = 0.325 

 puemis = puemis/exp((1.50e03-ctmax)/3.00e02) 

 uuemis = puemis 

c standard error cut off at impossible values follows 

 if (puemis.gt.(1.00-emissv)) puemis = 1.00-emissv 

 if (uuemis.gt.emissv) uuemis = emissv 

120 continue 

 return 

 end 

FRAPTRAN 1.4: ZOEMIS

*deck zoemis 

 subroutine zoemis (ctmax, zroxid, emissv) 

 implicit real (a-h,o-z) 

c 

c This is a MATPRO-11, Rev. 2 routine modified by PNNL for use in FRAPT 

c 

c zoemis calculates the emissivity of the cladding surface 

c as a function of maximum cladding temperature and 

c oxide thickness. expected standard error of the 

c calculated emissivity compared to in-reactor 

c data is also returned. 

c 

c emissv = output cladding surface emissivity (unitless) 

c puemis = output positive standard error expected in emissv 

c when compared to in-reactor data (not currently returned) 

c uuemis = output negative standard error expected in emissv 
A.25 

c when compared to in-reactor data (not currently returned) 

c 

c ctmax = input maximum cladding temperature (k) 

c zroxid = input oxide layer thickness (m) 

c 

c the equations used in this subroutine are based on data from 

c (1) aec fuels and materials development program progress 

c report no. 76, usaec report gemp - 1008 (1968). section by 

c e. f. juenke and s. j. sjodahl, p 239 

c (2) t. b. burgoyne and a. garlick, paper presented at 

c specialists meeting on the behaviour of water reactor 

c fuel elements under accident conditions, spatind norway 

c (september 1976) 

c (3) e. v. murphy and f. havelock, emissivity of zirconium 

c alloys in air in the temperature range 100 - 400 c 

c j. nuc. mat., 60 (1976) pp 167-176 

c 

c zoemis coded by r. l. miller sept 1974 

c modified by d. l. hagrman october 1976 

c Modified by PNNL, January 1997, to clean up coding, delete sensiti 

c uncertainty coding, and correct error in PUEMIS for T > 1500K 

c 

c model for temperatures below 1500K: 

 emissv = 3.25d-01 + 1.246d05*zroxid 

 if (zroxid .ge. 3.88d-06) emissv = 8.08642d-01 - 5.0d+01*zroxid 

 puemis = 0.1d0 

 uuemis = 0.1d0 

 if(ctmax .le. 1500.d0) go to 50 

c 

c modification for maximum temperatures above 1500K: 

 emissv = emissv * exp((1.50d03 - ctmax)/3.00d02) 

 if(emissv .lt. 0.325d0) emissv = 0.325d0 

 puemis = puemis*exp((1.50d03 - ctmax)/3.00d02) 

 uuemis = puemis 

c 

c standard error cut off at impossible values follows 

 if(puemis .gt. (1.00d0-emissv)) puemis = 1.00d0 - emissv 

 if(uuemis .gt. emissv) uuemis = emissv 

c 

 50 continue 

c 

 return 

c 

 end 

A.2.5 Cladding Thermal Expansion (CTHEXP) 

FRAPCON-3.4: CTHEXP

*deck cthexp 

 subroutine cthexp (ctemp,cathex,cdthex) 

c 

c cthexp calculates axial and diametral thermal expansion of 

c zircaloy lwr cladding. 

c cathex = output axial thermal expansion of zircaloy (m/m) 

c cdthex = output diametral thermal expansion of zircaloy (m/m) 

c ctemp = input cladding meshpoint temperature (k) 

c data for the correlation from room temperature to 1273k are 

c taken from r.l. mehan and f.w. wiesinger, "mechanical properties 

c of zircaloy-2", usaec report kapl-2110 (1961), d.b. scott, 

c "physical and mechanical properties of zircaloy-2 and -4", 

c wcap-3269-41 (1965), and j.j. kearns, "thermal expansion 
A.26 

c and preferred orientation in zircaloy", wapd-tm-472 (1965) 

c above 1273k the coefficient of thermal expansion used is the 

c constant value 9.7e-06/k recommended by b. lustman and f. kerze, 

c "the metallurgy of zirconium", mc graw-hill book company, new 

c york (1955) page 355. 

c between 1073 and 1273k (approximately the alpha-beta transition 

c range for zircaloy) cthexp uses the linear interpolation routine 

c polate to find the thermal expansion. 

c cthexp was originally coded by r.l. miller in nov. 1974. 

c last updated by g.a. reymann in june 1978. 

 dimension cathxd (44) 

 dimension dthexp (44) 

 data cathxd/ 0.0 , 293.15 , 

 + 3.52774e-03, 1073.15 , 3.53000e-03, 1083.15 , 

 + 3.50000e-03, 1093.15 , 3.46000e-03, 1103.15 , 

 + 3.41000e-03, 1113.15 , 3.33000e-03, 1123.15 , 

 + 3.21000e-03, 1133.15 , 3.07000e-03, 1143.15 , 

 + 2.80000e-03, 1153.15 , 2.50000e-03, 1163.15 , 

 + 2.00000e-03, 1173.15 , 1.50000e-03, 1183.15 , 

 + 1.30000e-03, 1193.15 , 1.16000e-03, 1203.15 , 

 + 1.13000e-03, 1213.15 , 1.10000e-03, 1223.15 , 

 + 1.11000e-03, 1233.15 , 1.13000e-03, 1243.15 , 

 + 1.20000e-03, 1253.15 , 1.30000e-03, 1263.15 , 

 + 1.40000e-03, 1273.15 / 

 data npcatx, iu / 22,1 / 

 data dthexp/ 4.3480e-04, 373.15, 

 + 5.1395e-03, 1073.15, 5.2200e-03, 1083.15, 5.2500e-03, 1093.15, 

 + 5.2800e-03, 1103.15, 5.2800e-03, 1113.15, 5.2400e-03, 1123.15, 

 + 5.2200e-03, 1133.15, 5.1500e-03, 1143.15, 5.0800e-03, 1153.15, 

 + 4.9000e-03, 1163.15, 4.7000e-03, 1173.15, 4.4500e-03, 1183.15, 

 + 4.1000e-03, 1193.15, 3.5000e-03, 1203.15, 3.1300e-03, 1213.15, 

 + 2.9700e-03, 1223.15, 2.9200e-03, 1233.15, 2.8700e-03, 1243.15, 

 + 2.8600e-03, 1253.15, 2.8800e-03, 1263.15, 2.9000e-03, 1273.15/ 

 data npdexp / 22 / 

 if (ctemp.le.1073.15) go to 100 

 if (ctemp.ge.1273.15) go to 110 

 cathex = polate(cathxd,ctemp,npcatx,iu) 

 cdthex = polate(dthexp,ctemp,npdexp,iu) 

 go to 120 

100 cathex = -2.5060e-05+(ctemp-273.15)*4.4410e-06 

 cdthex = -2.3730e-04+(ctemp-273.15)*6.7210e-06 

 go to 120 

110 cathex = -8.300e-03+(ctemp-273.15)*9.70e-06 

 cdthex = -6.800e-03+(ctemp-273.15)*9.70e-06 

120 continue 

 return 

 end 

FRAPTRAN 1.4: CTHEXP

*deck cthexp 

 subroutine cthexp (ctemp, cathex, cdthex, icm) 

 implicit real (a-h,o-z) 

c 

c cthexp calculates axial and diametral thermal expansion of 

c zircaloy lwr cladding. 

c 

c cathex = output axial thermal expansion of zircaloy (m/m) 

c cdthex = output diametral thermal expansion of zircaloy (m/m) 

c 

c ctemp = input cladding meshpoint temperature (k) 

c 
A.27 

c data for the correlation from room temperature to 1273k are 

c taken from r.l. mehan and f.w. wiesinger, "mechanical properties 

c of zircaloy-2", usaec report kapl-2110 (1961), d.b. scott, 

c "physical and mechanical properties of zircaloy-2 and -4", 

c wcap-3269-41 (1965), and j.j. kearns, "thermal expansion 

c and preferred orientation in zircaloy", wapd-tm-472 (1965) 

c 

c above 1273k the coefficient of thermal expansion used is the 

c constant value 9.7e-06/k recommended by b. lustman and f. kerze, 

c "the metallurgy of zirconium", mc graw-hill book company, new 

c york (1955) page 355. 

c 

c between 1073 and 1273k (approximately the alpha-beta transition 

c range for zircaloy) cthexp uses the linear interpolation routine 

c polate to find the thermal expansion. 

c 

c cthexp was originally coded by r.l. miller in nov. 1974. 

c last updated by g.a. reymann in june 1978. 

c Licensing analysis coding deleted 9/23/97 by ME Cunningham, PNNL 

c 

c generic Zry 

 dimension cathxd (44) 

 dimension dthexp (44) 

 data cathxd/ 0.0d0,293.15d0, 3.52774d-03,1073.15d0, 

 & 3.5300d-03,1083.15d0, 3.5000d-03,1093.15d0, 3.460d-03,1103.15d0, 

 & 3.4100d-03,1113.15d0, 3.3300d-03,1123.15d0, 3.210d-03,1133.15d0, 

 & 3.0700d-03,1143.15d0, 2.8000d-03,1153.15d0, 2.500d-03,1163.15d0, 

 & 2.0000d-03,1173.15d0, 1.5000d-03,1183.15d0, 1.300d-03,1193.15d0, 

 & 1.1600d-03,1203.15d0, 1.1300d-03,1213.15d0, 1.100d-03,1223.15d0, 

 & 1.1100d-03,1233.15d0, 1.1300d-03,1243.15d0, 1.200d-03,1253.15d0, 

 & 1.3000d-03,1263.15d0, 1.4000d-03,1273.15d0/ 

 data npcatx /22/ 

 data dthexp/ 4.3480d-04,373.15d0, 

 & 5.1395d-03,1073.15d0, 5.2200d-03,1083.15d0, 5.2500d-03,1093.15d0, 

 & 5.2800d-03,1103.15d0, 5.2800d-03,1113.15d0, 5.2400d-03,1123.15d0, 

 & 5.2200d-03,1133.15d0, 5.1500d-03,1143.15d0, 5.0800d-03,1153.15d0, 

 & 4.9000d-03,1163.15d0, 4.7000d-03,1173.15d0, 4.4500d-03,1183.15d0, 

 & 4.1000d-03,1193.15d0, 3.5000d-03,1203.15d0, 3.1300d-03,1213.15d0, 

 & 2.9700d-03,1223.15d0, 2.9200d-03,1233.15d0, 2.8700d-03,1243.15d0, 

 & 2.8600d-03,1253.15d0, 2.8800d-03,1263.15d0, 2.9000d-03,1273.15d0/ 

 data npdexp, iu /22, 1/ 

c 

 go to (400,400,400,400,400,600,400,600,400,400,400)icm 

c 

 400 continue 

c generic Zircaloy thermal expansion 

 if(ctemp .le. 1073.15d0) go to 401 

 if(ctemp .ge. 1273.15d0) go to 402 

c 

 cathex = polate(cathxd,ctemp,npcatx,iu) 

 cdthex = polate(dthexp,ctemp,npdexp,iu) 

 return 

c 

 401 continue 

 cathex = -2.5060d-05 + (ctemp-273.15d0)*4.4410d-06 

 cdthex = -2.3730d-04 + (ctemp-273.15d0)*6.7210d-06 

 return 

c 

 402 continue 

 cathex = -8.300d-03 + (ctemp-273.15d0)*9.70d-06 

 cdthex = -6.800d-03 + (ctemp-273.15d0)*9.70d-06 

 return 

c 
A.28 

 600 continue 

c RRC-KI thermal expansion 

c 

 if(ctemp .ge. 2133.d0) then 

 cathex=1.0582459d-2 

 cdthex=1.3133600d-2 

c 

 else if(ctemp .gt.1153.d0) then 

 cathex=1.076459d-3+9.7d-6*(ctemp-1153.d0) 

 cdthex=3.627600d-3+9.7d-6*(ctemp-1153.d0) 

c 

 else if(ctemp .gt. 883.d0)then 

 cathex=3.0465577d-3+2.312d-8*(ctemp-883.d0) 

 & -7.358d-8*(ctemp-883.d0)**2+1.7211d-10*(ctemp-883.d0)**3 

 cdthex=5.5977000d-3+2.312d-8*(ctemp-883.d0) 

 & -7.358d-8*(ctemp-883.d0)**2+1.7211d-10*(ctemp-883.d0)**3 

c 

 else if(ctemp .gt. 573.d0) then 

 cathex=0.13725577d-2+5.4d-6*(ctemp-573.d0) 

 cdthex=0.3336985d-8*ctemp**2+5.65390d-6*ctemp-0.199649865d-2 

c 

 else 

 cathex=0.1338985d-8*ctemp**2+3.85875d-6*ctemp-0.127813365d-2 

 cdthex=0.3336985d-8*ctemp**2+5.65390d-6*ctemp-0.199649865d-2 

 endif 

c 

 return 

 end 

A.2.6 Cladding Elastic Modulus (CELMOD) and Shear Modulus (CSHEAR) 

FRAPCON-3.4: CELMOD

*deck celmod 

 function celmod (ctemp,fnck,cwkf,deloxy) 

c 

c celmod calculates cladding young's modulus as a function 

c of temperature, fast neutron fluence, cold work, and average 

c oxygen concentration. grain orientation is assumed random. 

c the expected standard error is computed but not returned 

c (ucelmd). 

c celmod = output young's modulus for zircaloy 2 and 4 

c with random texture (pa) 

c ucelmd = output standard error expected in celmod (pa) 

c (not currently returned) 

c ctemp = input cladding temperature (k) 

c fnck = input effective fast fluence (neutrons/(m**2)) 

c cwkf = input effective cold work (unitless ratio of areas) 

c deloxy = input average oxygen concentration excluding 

c oxide layer - average oxygen concentration of 

c as-received cladding (kg oxygen/kg zircaloy) 

c the young's modulus calculated by this function is based on 

c data from the following references 

c (1) l. r. bunnell, g. b. mellinger and j. l. bates, high 

c temperature properties of zircaloy - oxygen alloys, epri 

c report np - 524 (1977). 

c (2) e. s. fisher and c. j. renken, physical review 135 2a 

c (20 july 1964) pp a482 - 494. 

c (3) p. e. armstrong and h. l. brown, transactions of the 

c metallurgical society of aime 230 (august 1964) 

c pp 962 - 966 

c (4) a. padel and a. groff, journal of nuclear materials 59 
A.29 

c (1976) pp 325-326 

c (5) w. r. smalley, saxton core ii fuel performance evaluation. 

c part i: materials, wcap-3385-56 (september 1971) 

c the expected standard error was derived using additional 

c data from the following references 

c (1) c. c. busby and c.r. woods (eds.) "properties of zircaloy-4 

c tubing", usaec report wapd-tm-585 (december 1966) p 65 

c (2) z. spasic, m. pavlovic and g. simis, conference on 

c the use of zirconium alloys in nuclear reactors, marianske 

c lanze, czech. conf-681086 (1968) pp 277 - 284 

c (3) r. l. mehan, modulus of elasticity of zircaloy-2 between 

c room temperature and 1000 f, kapl-m-rlm-16 (july 1958) 

c (4) d. o. northwood, i. m london, and l. e. bahen, journal 

c of nuclear materials 55 (1975) pp 299-310 

c (5) f. r. shober, j. a. van echo, l. l. marsh jr. and 

c j. r. keeler, the mechanical properties of zirconium and 

c zircaloy-2, bmi-1168 (1957) 

c celmod was coded by r. l. miller in march 1974. 

c updated and corrected by b. w. burnham in october 1975 

c modified by b. w. burnham october 1977 

c modified by d. l. hagrman december 1977 

 include 'lacmdl.h' 

 data on / 1 /, 

 + off / 2 /, 

 + locidx / 5 / 

 if (emflag(locidx).eq.on) go to 130 

c best estimate model young's modulus 

 c1 = (1.16e+11+ctemp*1.037e+08)*5.7015 

 c2 = 1.0 

 if (fnck.gt.1.0e+22) c2 = 0.88*(1.0-exp(-fnck/1.0e+25))+exp(-fnck/ 

 +1.0e+25) 

 c3 = -2.6e+10 

 celmod = (1.088e+11-5.475e+07*ctemp+c1*deloxy+c3*cwkf)/c2 

 if (ctemp.lt.1090.) go to 140 

c calculate a to a + b and a + b to b boundaries 

 wfox = deloxy+0.0012 

 taab = 1094.+wfox*(-1.289e+03+wfox*7.914e+05) 

 if (wfox.lt.0.025) go to 100 

 taab = 1556.4+3.8281e+04*(wfox-0.025) 

100 tabb = 392.46*((100*deloxy+0.1242807)**2+3.1417) 

 if (deloxy.lt.4.7308937e-03) go to 110 

 tabb = (100*deloxy+0.12)*491.157+1081.7413 

110 continue 

 if (ctemp.lt.taab) go to 140 

 if (ctemp.gt.tabb) go to 120 

 amodl = (1.088e+11-5.475e+07*taab+c1*deloxy+c3*cwkf)/c2 

 amodr = 9.21e+10-tabb*4.05e+07 

 celmod = amodl+(ctemp-taab)*(amodr-amodl)/(tabb-taab) 

 go to 140 

120 celmod = 9.21e+10-ctemp*4.05e+07 

 go to 140 

c em model young's modulus 

130 celmod = emclem(ctemp) 

140 continue 

 if (celmod.lt.1.0) celmod = 1.0 

 ucelmd = 6.4e+09 

 return 

 end 

FRAPTRAN1.4: CELMOD

*deck celmod 
A.30 

c 

 function celmod (ctemp, fnck, cwkf, deloxy, icm) 

c 

 implicit real (a-h,o-z) 

c 

c This is a MATPRO-11, Rev. 2 routine modified by PNNL for use in FRAPT 

c 

c celmod calculates cladding young's modulus as a function 

c of temperature, fast neutron fluence, cold work, and average 

c oxygen concentration. grain orientation is assumed random. 

c celmod = output young's modulus for zircaloy 2 and 4 

c with random texture (pa) 

c ctemp = input cladding temperature (k) 

c fnck = input effective fast fluence (neutrons/(m**2)) 

c cwkf = input effective cold work (unitless ratio of areas) 

c deloxy = input average oxygen concentration excluding 

c oxide layer - average oxygen concentration of 

c as-received cladding (kg oxygen/kg zircaloy) 

c 

c the young's modulus calculated by this function is based on 

c data from the following references 

c (1) l. r. bunnell, g. b. mellinger and j. l. bates, high 

c temperature properties of zircaloy - oxygen alloys, epri 

c report np - 524 (1977). 

c (2) e. s. fisher and c. j. renken, physical review 135 2a 

c (20 july 1964) pp a482 - 494. 

c (3) p. e. armstrong and h. l. brown, transactions of the 

c metallurgical society of aime 230 (august 1964) 

c pp 962 - 966 

c (4) a. padel and a. groff, journal of nuclear materials 59 

c (1976) pp 325-326 

c (5) w. r. smalley, saxton core ii fuel performance evaluation. 

c part i: materials, wcap-3385-56 (september 1971) 

c 

c the expected standard error was derived using additional 

c data from the following references 

c (1) c. c. busby and c.r. woods (eds.) "properties of zircaloy-4 

c tubing", usaec report wapd-tm-585 (december 1966) p 65 

c (2) z. spasic, m. pavlovic and g. simis, conference on 

c the use of zirconium alloys in nuclear reactors, marianske 

c lanze, czech. conf-681086 (1968) pp 277 - 284 

c (3) r. l. mehan, modulus of elasticity of zircaloy-2 between 

c room temperature and 1000 f, kapl-m-rlm-16 (july 1958) 

c (4) d. o. northwood, i. m london, and l. e. bahen, journal 

c of nuclear materials 55 (1975) pp 299-310 

c (5) f. r. shober, j. a. van echo, l. l. marsh jr. and 

c j. r. keeler, the mechanical properties of zirconium and 

c zircaloy-2, bmi-1168 (1957) 

c 

c celmod was coded by r. l. miller in march 1974. 

c updated and corrected by b. w. burnham in october 1975 

c modified by b. w. burnham october 1977 

c modified by d. l. hagrman december 1977 

c Modified by PNNL, January 1997, to clean up coding, delete licensi 

c analysis and sensitivity analysis coding 

c 

 go to (400,400,400,400,400,600,400,600,400,400,400)icm 

c 

 400 continue 

c generic Zry model 

c best estimate model young's modulus 

 c1 = (1.16d+11 + ctemp * 1.037d+08) * 5.7015d0 

 c2 = 1.0d0 
A.31 

 if(fnck .gt. 1.0d+22) c2 = 0.88d0 * (1.0d0 - exp(-fnck/1.0d+25)) 

 & + exp(-fnck/1.0d+25) 

 c3 = -2.6d+10 

 celmod = (1.088d+11 - 5.475d+07 * ctemp + c1 * deloxy + c3 * cwkf) 

 & /c2 

 if(ctemp .lt. 1090.d0) go to 401 

c 

c calculate a to a + b and a + b to b boundaries 

 wfox = deloxy + 0.0012d0 

 taab = 1094.d0 + wfox * (-1.289d+03 + wfox * 7.914d+05) 

 if (wfox .lt. 0.025d0) go to 415 

 taab = 1556.4d0 + 3.8281d+04 * (wfox - 0.025d0) 

 415 continue 

 tabb = 392.46d0 * ((100 * deloxy + 0.1242807d0)**2 + 3.1417d0) 

 if (deloxy .lt. 4.7308937d-03) go to 425 

 tabb = (100 * deloxy + 0.12d0) * 491.157d0 + 1081.7413d0 

 425 continue 

c 

 if (ctemp .lt. taab) go to 401 

 if (ctemp .gt. tabb) go to 435 

 amodl = (1.088d+11 - 5.475d+07 * taab + c1 * deloxy + c3 * cwkf) 

 & /c2 

 amodr = 9.21d+10 - tabb * 4.05d+07 

 celmod = amodl + (ctemp - taab) * (amodr - amodl)/(tabb - taab) 

 go to 401 

c 

 435 continue 

 celmod = 9.21d+10 - ctemp * 4.05d+07 

c 

 401 continue 

 if(celmod .lt. 1.0d0) celmod = 1.0d0 

 return 

c 

 600 continue 

c Zr-1%Nb model from RRC-KI 

c 

c This code block calculates cladding Young's modulus as 

c a function of temperature ; data is derived from 

c Volkov B.Yu. et.al."Material property library for 

c H1-alloy cladding", Preprint IAE-4941/11, 1989. 

c Fedorov et.al."Influence of oxygen content and gamma 

c -irradiation on the elastic module and internal 

c friction pure zirconium" Metal Physics vol.32(3) 

c p.621, 1971 

c 

 if(ctemp .gt. 1073.0d0) go to 601 

 celmod = 1.121d11 - 6.438d07*ctemp 

 celmod = celmod + 3.021d12*deloxy 

 go to 602 

 601 continue 

 celmod = 9.129d10 - 4.5d07*ctemp 

 602 continue 

 if(celmod .lt. 1.0d0) celmod = 1.0d0 

 return 

c 

 end 

FRAPCON-3.4: CSHEAR

*deck cshear 

 function cshear (ctemp,fnck,cwkf,deloxy) 

c 

c cshear calculates the shear modulus of zircaloy-2 and -4 
A.32 

c as a function of temperature, fast neutron fluence, 

c cold work, and average oxygen concentration. grain 

c orientation is assumed random. 

c the expected standard error is computed but not returned 

c (ucsher) 

c cshear = output shear modulus for zircaloy-2 and -4 

c with random texture (pa) 

c ucsher = output standard error expected in cshear (pa) 

c (not currently returned) 

c ctemp = input cladding temperature (k) 

c fnck = input effective fast fluence (neutrons/(m**2)) 

c cwkf = input effective cold work (unitless ratio of areas) 

c deloxy = input average oxygen concentration excluding 

c oxide layer - average oxygen concentration of 

c as-received cladding (kg oxygen/kg zircaloy) 

c the shear modulus calculated by this function is based on 

c data from the following references 

c (1) l. r. bunnell, g. b. mellinger and j. l. bates, high 

c temperature properties of zircaloy - oxygen alloys, epri 

c report np - 524 (1977) 

c (2) e. s. fisher and c. j. renken, physical review 135 2a 

c (20 july 1964) pp a482 - 494. 

c (3) p. e. armstrong and h. l. brown, transactions of the 

c metallurgical society of aime 230 (august 1964) 

c pp 962 - 966 

c (4) a. padel and a. groff, journal of nuclear materials 59 

c (1976) pp 325 - 326. 

c (5) w. r. smalley, saxton core ii fuel performance evaluation. 

c part i: materials, wcap-3385-56 (september 1971) 

c cshear was coded by r. l. miller, june 1974 

c modified by d. l. hagrman december 1977 

 c1 = 7.07e+11-ctemp*2.315e+08 

 c2 = 1.0 

 if (fnck.gt.1.0e+22) c2 = 0.88*(1.0-exp(-fnck/1.0e+25))+exp(-fnck/ 

 +1.0e+25) 

 c3 = -.867e+10 

 cshear = (4.04e+10-2.168e+07*ctemp+c1*deloxy+c3*cwkf)/c2 

 if (ctemp.lt.1090.) go to 130 

c calculate a to a + b and a + b to b boundaries 

 wfox = deloxy+0.0012 

 taab = 1094.+wfox*(-1.289e+03+wfox*7.914e+05) 

 if (wfox.lt.0.025) go to 100 

 taab = 1556.4+3.8281e+04*(wfox-0.025) 

100 tabb = 392.46*((100*deloxy+0.1242807)**2+3.1417) 

 if (deloxy.lt.4.7308937e-03) go to 110 

 tabb = (100*deloxy+0.12)*491.157+1081.7413 

110 continue 

 if (ctemp.lt.taab) go to 130 

 if (ctemp.gt.tabb) go to 120 

 amodl = (4.04e+10-2.168e+07*taab+c1*deloxy+c3*cwkf)/c2 

 amodr = 3.49e+10-tabb*1.66e+07 

 cshear = amodl+(ctemp-taab)*(amodr-amodl)/(tabb-taab) 

 go to 130 

120 cshear = 3.49e+10-ctemp*1.66e+07 

130 continue 

 if (cshear.lt.1.0) cshear = 1.0 

 ucsher = 9.0e+09 

 return 

 end 

FRAPTRAN 1.4: CSHEAR
A.33 

*deck cshear 

c 

 function cshear (ctemp, fnck, cwkf, deloxy, icm) 

c 

 implicit real (a-h,o-z) 

c 

c This is a MATPRO-11, Rev. 1 routine modified by PNNL for use in FRAPT 

c 

c cshear calculates the shear modulus of zircaloy-2 and -4 

c as a function of temperature, fast neutron fluence, 

c cold work, and average oxygen concentration. grain 

c orientation is assumed random. 

c cshear = output shear modulus for zircaloy-2 and -4 

c with random texture (pa) 

c ctemp = input cladding temperature (k) 

c fnck = input effective fast fluence (neutrons/(m**2)) 

c cwkf = input effective cold work (unitless ratio of areas) 

c deloxy = input average oxygen concentration excluding 

c oxide layer - average oxygen concentration of 

c as-received cladding (kg oxygen/kg zircaloy) 

c 

c the shear modulus calculated by this function is based on 

c data from the following references 

c 

c (1) l. r. bunnell, g. b. mellinger and j. l. bates, high 

c temperature properties of zircaloy - oxygen alloys, epri 

c report np - 524 (1977) 

c (2) e. s. fisher and c. j. renken, physical review 135 2a 

c (20 july 1964) pp a482 - 494. 

c (3) p. e. armstrong and h. l. brown, transactions of the 

c metallurgical society of aime 230 (august 1964) 

c pp 962 - 966 

c (4) a. padel and a. groff, journal of nuclear materials 59 

c (1976) pp 325 - 326. 

c (5) w. r. smalley, saxton core ii fuel performance evaluation. 

c part i: materials, wcap-3385-56 (september 1971) 

c 

c cshear was coded by r. l. miller, june 1974 

c modified by d. l. hagrman december 1977 

c Modified by PNNL, January 1997, to clean up coding and delete 

c sensitivity uncertainty analysis coding 

c 

 go to (400,400,400,400,400,600,400,600,400,400,400)icm 

c 

 400 continue 

c generic Zry properties 

 c1 = 7.07d+11 - ctemp * 2.315d+08 

 c2 = 1.0d0 

 if(fnck .gt. 1.0d+22) c2 = 0.88d0 * (1.0d0 - exp(-fnck/1.0d+25)) 

 & + exp(-fnck/1.0d+25) 

 c3 = -0.867d+10 

 cshear = (4.04d+10 - 2.168d+07 * ctemp + c1 * deloxy + c3*cwkf)/c2 

 if(ctemp .lt. 1090.d0) go to 401 

c 

c calculate a to a + b and a + b to b boundaries 

 wfox = deloxy + 0.0012d0 

 taab = 1094.d0 + wfox * (-1.289d+03 + wfox * 7.914d+05) 

 if(wfox .lt. 0.025d0) go to 415 

 taab = 1556.4d0 + 3.8281d+04 * (wfox - 0.025d0) 

 415 continue 

 tabb = 392.46d0 * ((100 * deloxy + 0.1242807d0)**2 + 3.1417d0) 

 if(deloxy .lt. 4.7308937d-03) go to 425 

 tabb = (100 * deloxy + 0.12d0) * 491.157d0 + 1081.7413d0 
A.34 

 425 continue 

c 

 if(ctemp .lt. taab) go to 401 

 if(ctemp .gt. tabb) go to 435 

 amodl = (4.04d+10 - 2.168d+07 * taab + c1 * deloxy + c3*cwkf)/c2 

 amodr = 3.49d+10 - tabb * 1.66d+07 

 cshear = amodl + (ctemp - taab ) * (amodr - amodl )/(tabb - taab) 

 go to 401 

 435 continue 

 cshear = 3.49d+10 - ctemp * 1.66d+07 

 401 continue 

 if(cshear .lt. 1.0d0) cshear = 1.0d0 

 return 

c 

 600 continue 

c Zr-1%Nb from RRC-KI 

c 

c This code block calculates cladding young's modulus and 

c poisson's ratio as a function of temperature ; data is 

c derived from Volkov B.Yu. et.al."Material property library 

c for H1-alloy cladding", Preprint IAE-4941/11, 1989. 

c 

c elastic modulus 

 ey = celmod (ctemp, fnck, cwkf, deloxy, CladType) 

c 

c poisson's ratio: 

 xnu = 0.42628d0 - 5.556d-5*ctemp 

c 

c shear modulus 

 cshear = 0.5d0*ey/(1.d0+xnu) 

 if(cshear .lt. 1.0d0) cshear = 1.0d0 

 return 

c 

 end 

A.2.7 Cladding Axial Growth (CAGROW) 

FRAPCON-3.4: CAGROW

*deck cagrow 

 function cagrow (flux,fluenc,icm,dtime) 

 common /uncert/sigftc, sigftex, sigfgr 

 & ,sigswell ,sigcreep ,siggro ,sigcor ,sigh2 

c 

c cagrow calculates axial cladding growth strain increment as a functi 

c flux = input fast neutron flux ((neutrons/m**2)/s) 

c fluenc = input fast neutron fluence (neutrons/m**2) 

c dtime = input time increment at temperature (s) 

c icm = input cladding material index:4 - Zircaloy 4 

c 2 - Zircaloy 2 

c cagrow = output axial cladding growth straini increment (m/m) 

c cagrow coded by D. L. Hagrman, August 1975 

c revised April, 1995 by D. D. Lanning and K. J. Geelhood 

c to conform to EPRI model by D. G. Franklin, "Zircaloy-4 

c Cladding Deformation During Power Reactor Irradiation," 

c ASTM STP 754, 1982, pp.235-267. 

 f2=fluenc/10000. 

 f1=(fluenc-flux*dtime)/10000. 

 ax1=2.18e-21*f1**0.845 

 ax2=2.18e-21*f2**0.845 

 if(icm.eq.5) then 

 ax1=7.013e-21*f1**0.81787 
A.35 

 ax2=7.013e-21*f2**0.81787 

 elseif(icm.eq.6) then 

 ax1=9.7893e-25*f1**0.98239 

 ax2=9.7893e-25*f2**0.98239 

 endif 

 cagrow = ax2-ax1 

 if( icm .lt. 4.0 ) cagrow = cagrow * 0.5 

c add on uncertainty 

c PWR cladding 

 if(siggro.gt.0.0.and.icm.ge.4) 

 &cagrow = cagrow*(1.0+siggro*0.223) 

 if(siggro.lt.0.0.and.icm.ge.4) 

 &cagrow = cagrow/(1.0-siggro*0.223) 

c BWR cladding 

 if(siggro.gt.0.0.and.icm.lt.4) 

 &cagrow = cagrow*(1.0+siggro*0.203) 

 if(siggro.lt.0.0.and.icm.lt.4) 

 &cagrow = cagrow/(1.0-siggro*0.203) 

 return 

 end 

A.2.8 Cladding Creep Rate (CREPR) 

FRAPCON-3.4: CREPR

*deck crepr 

c 

 subroutine crepr(sig,edot,phi,fluence,cldwks,CladAveTemp 

 &,GapPress,CoolantPress,rci,rco, ProblemTime) 

c 

 common /cladtype/ icm 

 common /uncert/sigftc, sigftex, sigfgr 

 & ,sigswell ,sigcreep ,siggro ,sigcor ,sigh2 

c implicit real*8 (a-h,o-z) 

c ******************************************************************** 

c crepr is called from cladf. 

c crepr computes (transverse) cladding creep strain rate 

c as a function of effective stress,temperature, fast 

c neutron flux and accumulated total creep strain 

c ******************************************************************** 

c input arguments 

c ******************************************************************** 

c CreepStrain - total accumulated creep strain (dimensionless) 

c phi - flux (neutrons/m**2/sec) 

c ProblemTime - time to end of power step (hr) 

c CladAveTemp - cladding average temperature (F) 

c sig - cladding hoop stress per node (psi) 

c ******************************************************************** 

c output arguments 

c ******************************************************************** 

c edot - creep strain rate (1/hr) 

c ******************************************************************** 

c 

c NEW MODEL 

c New model uses a modified version of the Limback and Andersson equation 

c ASTM STP 1295 pp. 448-468 

c calculate stress in each direction 

 sr=(GapPress*rci**2-CoolantPress*rco**2 

 &+rci**2*rco**2*(CoolantPress-GapPress)/((rco+rci)/2.0)**2) 

 &/(rco**2-rci**2) 

 st=(GapPress*rci**2-CoolantPress*rco**2 

 &-rci**2*rco**2*(CoolantPress-GapPress)/((rco+rci)/2.0)**2) 
A.36 

 &/(rco**2-rci**2) 

 sl=(GapPress*rci**2-CoolantPress*rco**2)/(rco**2-rci**2) 

c sag is effective stress in MPa 

c 

c sag = sig*6894.7573/1.0E6 

 sag=0.707*sqrt((sl-st)**2+(st-sr)**2+(sr-sl)**2)*6894.7573/1.0E6 

 tcak = (CladAveTemp+459.67)/1.8 

 fluenceLimback = fluence/10000.0 

c 

c Thermal Creep Rate 

c 

 A = 1.08E9 

 if(cldwks.eq.0.0.or.icm.eq.5) A = 5.47E8 

 an = 2.0 

 if(cldwks.eq.0.0.or.icm.eq.5) an = 3.5 

 E = 1.148E5 - 59.9*tcak 

 Q = 201000.0 

 R = 8.314 

 ai = 650.0*(1.0-0.56*(1.0-exp(-1.4e-27*fluenceLimback**1.3))) 

 edottherm = A*(E/tcak)*(sinh(ai*abs(sag)/E))**an*exp(-Q/(R*tcak)) 

c 

c Irradiation Creep Rate 

c 

 c0SR = 4.0985E-24 

 c0RXA = 1.87473E-24 

 c0 = C0SR 

 if(cldwks.eq.0.0) c0 = c0RXA 

 c1 = 0.85 

 c2 = 1.0 

 tmod=tcak 

 if(tcak.lt.570.0) tmod=570.0 

 if(tcak.gt.625.0) tmod=625.0 

 ctemp=-7.0237+0.0136*tmod 

 if(cldwks.eq.0.0) ctemp=-3.18562+0.00699132*tmod 

 edotirr = c0*phi**c1*(abs(sag))**c2*ctemp 

c 

c Primary Creep Rate 

c 

 if(ProblemTime.ge.0.1) then 

 esat=0.0216*(edottherm+edotirr)**0.109* 

 &(2.0-tanh(3.55e4*(edottherm+edotirr)))**-2.05 

 edotprimary= 

 &0.5*esat*52.0*(edottherm+edotirr)**0.5/((ProblemTime)**0.5)* 

 &exp(-52.0*(edottherm+edotirr)**0.5*(ProblemTime)**0.5) 

 edottherm=edottherm+edotprimary 

 endif 

c 

 edot = edottherm + edotirr 

 edot = edot*1.1547 

 if(icm.eq.6) edot = edot*0.8 

100 continue 

c add on uncertainty 

c SRA cladding 

 if(sigcreep.gt.0.0.and.cldwks.ne.0.0) 

 &edot = edot*(1.0+sigcreep*0.145) 

 if(sigcreep.lt.0.0.and.cldwks.ne.0.0) 

 &edot = edot/(1.0-sigcreep*0.145) 

c RXA cladding 

 if(sigcreep.gt.0.0.and.cldwks.eq.0.0) 

 &edot = edot*(1.0+sigcreep*0.216) 

 if(sigcreep.lt.0.0.and.cldwks.eq.0.0) 

 &edot = edot/(1.0-sigcreep*0.216) 

 return 
A.37 

 end 

A.2.9 Cladding Meyer Hardness (CMHARD) 

FRAPCON-3.4: CMHARD

*deck cmhard 

 function cmhard (ctemp) 

c 

c the routine cmhard calculates meyer hardness as a function of 

c cladding temperature. 

c cmhard = output meyer hardness of zircaloy cladding (n/m**2) 

c ctemp = input cladding temperature (k) 

c the equations used in this function are based on data from 

c (1) a. m. ross and r. l. stoute, heat transfer coefficient 

c between uo2 and zircaloy - 2, aecl - 1552 (june 1962) 

c (2) i. d. peggs and d. p. godin, the yield strength - hot 

c hardness relationship of zircaloy - 4, journal of nuclear 

c materials 57 pp 246 - 248 (1975) 

c cmhard was coded by v.f.baston in may 1974. 

c modified by m. a. morgan june 1978 

 t = ctemp 

 cmhard = exp(2.6034e01-2.6394e-02*t+4.3502e-05*t**2-2.5621e-08*t** 

 +3) 

 return 

 end 

FRAPTRAN 1.4: CMHARD

*deck cmhard 

 function cmhard (ctemp, icm) 

 implicit real (a-h,o-z) 

c 

c This is a MATPRO-11, Rev. 2 routine modified by PNNL for use in FRAPT 

c 

c the routine cmhard calculates meyer hardness as a function of 

c cladding temperature. 

c 

c cmhard = output meyer hardness of zircaloy cladding (n/m**2) 

c 

c ctemp = input cladding temperature (k) 

c 

c the equations used in this function are based on data from 

c (1) a. m. ross and r. l. stoute, heat transfer coefficient 

c between uo2 and zircaloy - 2, aecl - 1552 (june 1962) 

c (2) i. d. peggs and d. p. godin, the yield strength - hot 

c hardness relationship of zircaloy - 4, journal of nuclear 

c materials 57 pp 246 - 248 (1975) 

c 

c cmhard was coded by v.f.baston in may 1974. 

c modified by m. a. morgan june 1978 

c Modified by PNNL, January 1997, to clean up coding, delete sensiti 

c uncertainty coding 

c 

 go to (400,400,400,400,400,600,400,600,400,400,400)icm 

c 

 400 continue 

c generic Zry property 

 cmhard = exp(2.6034d+01 - 2.6394d-02*ctemp + 4.3502d-05*ctemp**2 

 & - 2.5621d-08*ctemp**3) 

 if (cmhard .lt. 1.94d+08) cmhard = 1.94d+08 
A.38 

 return 

c 

 600 continue 

c Zr-1%Nb from RRC-KI 

 if(ctemp .lt. 800.d0) cmhard=1.d6*(2172.1d0 - 10.7055d0*ctemp 

 & + 0.02765d0*ctemp**2 - 3.278d-5*ctemp**3 + 1.423d-8*ctemp**4) 

 if(ctemp .ge. 800.d0) cmhard=exp(26.034d0 - 2.6394d-2*ctemp 

 & + 4.3502d-5*ctemp**2 - 2.5621d-8*ctemp**3) 

 if (cmhard .lt. 1.0d+05) cmhard = 1.0d+05 

c 

 end 

A.3 Gas Material Properties 

A.3.1 Gas Thermal Conductivity (GTHCON) 

FRAPCON-3.4

*deck gthcon 

 function gthcon (gmix,gtemp,gpres,gpthk) 

c 

c gthcon calculates gas thermal conductivity as a function of 

c temperature and gas fraction for seven gases: 

c gthcon = output gas thermal conductivity (w/m-k). 

c gmix(i) = input mole fractions of the gas mixture 

c the seven elements of gmix must sum to 1.0 

c constituent gas number key 

c 1 helium 

c 2 argon 

c 3 krypton 

c 4 xenon 

c 5 hydrogen 

c 6 nitrogen 

c 7 water vapor 

c gtemp = input gas temperature (k) 

c gpres = input gas pressure (pa) 

c used for knudsen domain correction and for steam 

c not used if 0. 

c gpthk = input effective gap thickness for knudsen domain (m) 

c (maximum of gap dimension or surface roughness) 

c not used if 0. 

c formula for gas mixtures is from r.s.brokaw, report nasa tr r-81 

c (1960). conductivity of rare gases is based on j.m.gandhi and 

c s.c.saxena, jour. chem. and eng. data, vol.13, no.3 (1968) 

c also: wisconsin electric power co., docket no. 50-301 (jan 1973) 

c the accommodation factor is from r.a.dean, cvna-127 (1962) 

c steam equation is from meyer, et.al., "thermodynamic and 

c transport properties of steam", the american society of 

c mechanical engineers (1967) 

c gthcon coded by r.c.young march 1975 

c adapted from routine cmix by p.e.macdonald 

 include 'lacmdl.h' 

 dimension gmix(7), a(7), c(7), r(7), aa(6), bb(6) 

 data a /4.003 , 39.944 ,83.80 ,131.30 , 2.016 ,28.8 , 

 + 18.016 / 

c r = sqrt(a) 

 data r /2.00075 , 6.3201 , 9.1542 ,11.4586 , 1.41986 , 

 + 5.3666 , 4.2445 / 

c 

c The "aa" and "bb" parameters are from MATPRO-11 Rev. 1, 1981 

c "aa" and "bb" parameters updated by W.G. Luscher 10-17-07 
A.39 

c 

 data aa/2.531e-3,4.092e-4,1.966e-4,9.825e-5,1.349e-3,2.984e-4/ 

 data bb/0.7146, 0.6748, 0.7006, 0.7334, 0.8408, 0.7799/ 

c 

c 

 data eps/1.e-9 / 

 data on / 1 /, 

 + off / 2 /, 

 + locidx / 16 / 

 if (emflag(locidx).eq.on) gthcon = emgton(gmix,gtemp,gpres,gpthk) 

 if (emflag(locidx).eq.on) go to 170 

 if (gtemp.gt.5555.) go to 100 

 if (gtemp.gt.0.0000) go to 110 

100 write (6,180) gtemp,gmix 

 stop 

110 l = 0 

c 

 do 120 i=1,6 

 if (gmix(i).le.eps) go to 120 

 l = i 

 c(i) = aa(i)*gtemp**bb(i) 

120 continue 

 if (gmix(7).le.eps) go to 140 

 l = 7 

c steam 

c c(7) defination updated by W.G. Luscher 10-17-07 

 tc = gtemp - 273.15d0 

 if (gtemp.le.973.15d0) 

 & c(7) = gpres/gtemp*(-2.8516d-08 + 9.424d-10*gtemp 

 & -6.005d-14*gtemp**2) + 1.009d0*gpres**2/gtemp**2 

 & /(tc)**4.2d0 + (17.6d-4 + 5.87d-5*(tc) 

 & + 1.08d-7*(tc)**2 - 4.51d-11*(tc)**3) 

 if (gtemp.gt.973.15d0) 

 & c(7) = 4.44d-6*gtemp**1.45d0 + 9.45d0-05*(2.1668d-09*gpres 

 & /gtemp)**1.3d0 

140 if (l.le.0) go to 100 

 gthcon = 0. 

 do 160 i=1,l 

 if (gmix(i).le.eps) go to 160 

 sum = 0. 

 do 150 j=1,l 

 if (j.eq.i) go to 150 

 if (gmix(j).le.eps) go to 150 

 rc = c(i)/c(j) 

 ra = a(i)/a(j) 

 fij = 1.+2.41*(ra-1.)*(ra-.142)/(1.+ra)**2 

 gij = (1.+sqrt(rc*r(i)/r(j)))**2/sqrt(8.*(1.+ra)) 

 sum = sum+fij*gij*gmix(j) 

150 continue 

 gthcon = gthcon+c(i)*gmix(i)/(gmix(i)+sum) 

160 continue 

170 continue 

 return 

180 format ('1gthcon input bad. gtemp=',1pe12.4,' k',/'0fractions',( 

 +10e12.4)) 

 end 

FRAPTRAN 1.4

*deck gthcon 

 function gthcon (gmix,gtemp,gpres,gpthk) 

 implicit real (a-h, o-z) 
A.40 

c 

c gthcon calculates gas thermal conductivity as a function of 

c temperature and gas fraction for seven gases: 

c 

c gthcon = output gas thermal conductivity (w/m-K). 

c gmix(i) = input mole fractions of the gas mixture 

c the seven elements of gmix must sum to 1.0 

c constituent gas number key 

c 1 helium 

c 2 argon 

c 3 krypton 

c 4 xenon 

c 5 hydrogen 

c 6 nitrogen 

c 7 water vapor 

c 

c gtemp = input gas temperature (K) 

c gpres = input gas pressure (pa) 

c used for knudsen domain correction and for steam 

c not used if 0. 

c gpthk = input effective gap thickness for knudsen domain (m) 

c (maximum of gap dimension or surface roughness) 

c not used if 0. 

c 

c formula for gas mixtures is from r.s.brokaw, report nasa tr r-81 

c (1960). conductivity of rare gases is based on j.m.gandhi and 

c s.c.saxena, jour. chem. and eng. data, vol.13, no.3 (1968) 

c also: wisconsin electric power co., docket no. 50-301 (jan 1973) 

c the accommodation factor is from r.a.dean, cvna-127 (1962) 

c steam equation is from meyer, et.al., "thermodynamic and 

c transport properties of steam", the american society of 

c mechanical engineers (1967) 

c 

c gthcon coded by r.c.young march 1975 

c adapted from routine cmix by p.e.macdonald 

c 

 dimension gmix(7), a(7), c(7), r(7), aa(6), bb(6) 

 data a /4.003d0, 39.944d0, 83.80d0, 131.30d0, 2.016d0, 

 & 28.02d0, 18.02d0/ 

c r = sqrt(a) 

 data r /2.00075d0, 6.3201d0, 9.1542d0, 11.4586d0, 1.41986d0, 

 & 5.3666d0, 4.2445d0 / 

c 

c The "aa" and "bb" parameters are from MATPRO-11 Rev. 2, 8/1981 

c "aa" and "bb" parameters updated by W.G. Luscher 10-17-07 

c 

 data aa/2.531d-3,4.092d-4,1.966d-4,9.825d-5,1.349d-3,2.984d-4/ 

 data bb/0.7146d0,0.6748d0,0.7006d0,0.7334d0,0.8408d0,0.7799d0/ 

 data eps/1.d-9 / 

c 

 if (gtemp.gt.5555.d0) go to 100 

 if (gtemp.gt.0.0000d0) go to 110 

 100 continue 

 write (6,180) gtemp,gmix 

 stop 

c 

 110 continue 

 l = 0 

c 

 do 120 i=1,6 

 if (gmix(i).le.eps) go to 120 

 l = i 

 c(i) = aa(i)*gtemp**bb(i) 
A.41 

 120 continue 

c 

 if (gmix(7).le.eps) go to 140 

 l = 7 

c steam (from MATPRO-11, Rev. 2) 

c c(7) definition updated by W.G. Luscher 10-17-07 

 tc = gtemp -273.15d0 

 if (gtemp .le. 973.15d0) 

 & c(7) = gpres/gtemp*(-2.8516d-08 + 9.424d-10*gtemp 

 & - 6.005d-14*gtemp**2) + 1.009d0*gpres**2/gtemp**2 

 & /(gtemp-273.15d0)**4.2d0 + (17.6d-4 + 5.87d-5*(tc) 

 & + 1.08d-7*(tc)**2 - 4.51d-11*(tc)**3) 

 if (gtemp .gt. 973.15d0) 

 & c(7) = 4.44d-06*gtemp**1.45d0 + 9.45d-05*(2.1668d-09*gpres 

 & /gtemp)**1.3d0 

c c(7) = 17.6e-3+tc*(5.87e-5+tc*(1.04e-7-4.51e-11*tc)) 

 140 continue 

 if (l.le.0) go to 100 

 gthcon = 0.d0 

c 

 do 160 i=1,l 

 if (gmix(i).le.eps) go to 160 

 sum = 0.d0 

 do 150 j=1,l 

 if (j.eq.i) go to 150 

 if (gmix(j).le.eps) go to 150 

 rc = c(i)/c(j) 

 ra = a(i)/a(j) 

 fij = 1.d0+2.41d0*(ra-1.d0)*(ra-.142d0)/(1.d0+ra)**2 

 gij = (1.d0+sqrt(rc*r(i)/r(j)))**2/sqrt(8.d0*(1.d0+ra)) 

 sum = sum+fij*gij*gmix(j) 

 150 continue 

 gthcon = gthcon+c(i)*gmix(i)/(gmix(i)+sum) 

 160 continue 

c 

 return 

c 

 180 format ('gthcon input bad. gtemp= ',1pe12.4,' K',/' fractions', 

 & (10e12.4)) 

c 

 end