Material Property Correlations: 

Comparisons between 

FRAPCON-3.5, FRAPTRAN-1.5, 

and MATPRO

Office of Nuclear Regulatory Research

NUREG/CR-7024, Rev. 1

PNNL-19417, Rev. 1
NRC Reference Material 

As of November 1999, you may electronically access 

NUREG-series publications and other NRC records at 

NRC’s Public Electronic Reading Room at 

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

records include, to name a few, NUREG-series 

publications; Federal Register notices; applicant, 

licensee, and vendor documents and correspondence; 

NRC correspondence and internal memoranda; bulletins 

and information notices; inspection and investigative 

reports; licensee event reports; and Commission papers 

and their attachments. 

NRC publications in the NUREG series, NRC 

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

Federal Regulations may also be purchased from one 

of these two sources. 

1. The Superintendent of Documents

 U.S. Government Printing Office

Mail Stop SSOP

Washington, DC 20402–0001

Internet: bookstore.gpo.gov

Telephone: 202-512-1800

Fax: 202-512-2250 

2. The National Technical Information Service 

Springfield, VA 22161–0002

www.ntis.gov

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

A single copy of each NRC draft report for comment is 

available free, to the extent of supply, upon written 

request as follows:

Address: U.S. Nuclear Regulatory Commission

 Office of Administration 

Publications Branch

 Washington, DC 20555-0001

E-mail: DISTRIBUTION.RESOURCE@NRC.GOV

Facsimile: 301–415–2289 

Some publications in the NUREG series that are 

posted at NRC’s Web site address 

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

are updated periodically and may differ from the last 

printed version. Although references to material found on 

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

material available on the date cited may subsequently be 

removed from the site.

Non-NRC Reference Material 

Documents available from public and special technical 

libraries include all open literature items, such as books, 

journal articles, transactions, Federal Register notices, 

Federal and State legislation, and congressional reports. 

Such documents as theses, dissertations, foreign reports 

and translations, and non-NRC conference proceedings 

may be purchased from their sponsoring organization. 

Copies of industry codes and standards used in a 

substantive manner in the NRC regulatory process are 

maintained at—

The NRC Technical Library

Two White Flint North

11545 Rockville Pike

Rockville, MD 20852–2738 

These standards are available in the library for reference 

use by the public. Codes and standards are usually 

copyrighted and may be purchased from the originating 

organization or, if they are American National Standards, 

from—

American National Standards Institute

11 West 42nd Street

New York, NY 10036–8002

www.ansi.org

212–642–4900 

AVAILABILITY OF REFERENCE MATERIALS

IN NRC PUBLICATIONS

Legally binding regulatory requirements are stated only 

in laws; NRC regulations; licenses, including technical 

specifications; or orders, not in NUREG-series 

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

those of the NRC.

The NUREG series comprises (1) technical and 

administrative reports and books prepared by the staff 

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

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

XXXX), (3) reports resulting from international 

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

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

decisions and orders of the Commission and Atomic and 

Safety Licensing Boards and of Directors’ decisions

under Section 2.206 of NRC’s regulations (NUREG–

0750).

DISCLAIMER: This report was prepared as an account 

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

Government. Neither the U.S. Government nor any 

agency thereof, nor any employee, makes any warranty, 

expressed or implied, or assumes any legal liability or 

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

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

process disclosed in this publication, or represents that 

its use by such third party would not infringe privately 

owned rights.
Material Property Correlations: 

Comparisons between 

FRAPCON-3.5, FRAPTRAN-1.5, 

and MATPRO

Manuscript Completed: May 2014

Date Published: October 2014

Prepared by:

W. G. Luscher and K. J. Geelhood

Pacific Northwest National Laboratory

Richland, WA 99352

P. Raynaud, NRC Project Manager

NRC Job Code V6197

Office of Nuclear Regulatory Research

NUREG/CR-7024, Rev. 1

PNNL-19417, Rev. 1

Abstract

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

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

analyses. 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.5 and FRAPTRAN-1.5. 

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.5 and FRAPTRAN-1.5 correlations are provided.

 

iii

Foreword

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

to model steady-state and transient fuel behavior, respectively, to support regulatory analyses. 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.5 and FRAPTRAN-1.5, this report also provides a variety of comparisons to the material 

property correlations used in MATPRO.

v

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-27

2.5.3 References............................................................................................................... 2-27

2.6 Fuel Densification (FUDENS)........................................................................................ 2-27

2.6.1 Model Development and Comparisons................................................................... 2-28

2.6.2 Applicability and Uncertainty................................................................................. 2-31

2.6.3 References............................................................................................................... 2-31

2.7 Fuel Swelling (FSWELL)............................................................................................... 2-31

2.7.1 Model Development and Comparisons................................................................... 2-32

2.7.2 Applicability and Uncertainty................................................................................. 2-35

2.7.3 References............................................................................................................... 2-35

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-4

3.2.1 Model Development and Comparisons..................................................................... 3-4

3.2.2 Applicability and Uncertainty................................................................................... 3-7

3.2.3 References................................................................................................................. 3-7

3.3 Cladding Oxide Thermal Conductivity (ZOTCON)......................................................... 3-8

vii
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-8

4.1.3 References................................................................................................................. 4-8

Appendix A Subroutine Source Codes........................................................................................ A-1

viii
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 

unirradiated 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 unirradiated MOX 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-26

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-34

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-35

ix
Figure 3.2-1 CTHCON model and data used for development. .......................................... 3-6

Figure 3.2-2 CTHCON model and data collected since development................................. 3-6

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

temperature compared with the MATPRO and FRAPCON-3.5/

FRAPTRAN-1.5 (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.5/FRAPTRAN-1.5 

(FRAP) and MATPRO axial thermal expansion. Data from cladding 

tubes. ............................................................................................................. 3-17

Figure 3.5-2 Model-to-data comparison for FRAPCON-3.5/FRAPTRAN-1.5 

(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

x
List of Tables

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

performance codes, FRAPCON-3.5 and FRAPTRAN-1.5............................. 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-25

Table 3.0-1 Example heat treatments and coldworked conditions for different alloys. ..... 3-1

Table 3.1-1 Specific heat capacity database for Zircaloy-2, Zircaloy-4, ZIRLO, 

Optimized ZIRLO, and M5............................................................................. 3-2

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

Table 3.8-1 Parameters for FRAPCON-3.5 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

xi

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, to support regulatory 

analyses. 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.5 and FRAPTRAN-1.5. 

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.5 and FRAPTRAN-1.5 correlations are also provided.

This is one of three documents that describe the FRAPCON-3.5 and FRAPTRAN-1.5 fuel 

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

the FRAPCON-3.5 code description document (Geelhood and Luscher, 2014). The FRAPTRAN1.5 code structure and behavioral models are described in the FRAPTRAN-1.5 code description 

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

report (Geelhood, 2008).

xiii

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

xv
Xe xenon 

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

ZrO2 zirconium dioxide

μm micrometer(s)

xvi
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, to support regulatory

analyses. 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.5 and FRAPTRAN1.5. The source code for each material property correlation discussed is provided for FRAPCON3.5 and FRAPTRAN-1.5 (see Appendix A) as well as a range of applicability and an estimate of 

uncertainty, where possible.

1-1
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.5

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

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

description document (Geelhood and Luscher, 2014). The FRAPTRAN-1.5 code structure and 

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

al., 2014). 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.5 and FRAPTRAN-1.5. 

Model/Property FRAPCON-3.5 FRAPTRAN-1.5

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 – solid Material properties handbook NA

Fuel swelling – gaseous 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

1-2
Model/Property FRAPCON-3.5 FRAPTRAN-1.5

Cladding specific heat Material properties handbook Material properties handbook

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

(Knuutila, 2006)

VTT-R-11337-06

(Knuutila, 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-3
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. and W.G. Luscher. 2014. 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 Rev.1, PNNL-19418, Vol. 1 Rev.1, Pacific Northwest National 

Laboratory, Richland, WA. 

Geelhood, K.J., W.G. Luscher, and J.M. Cuta. 2014. FRAPTRAN 1.4: A Computer Code for the 

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

Rev.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.

1-4
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.5 and FRAPTRAN-1.5 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)

3 2 sldus(c) 2840 5.41395c 7.468390x10 c − = − +

2 2 liqdus(c) 2840 3.21860c 1.448518x10 c − = − +

ftmelt = 3113.15 − 5FBu /10000

fdelta = liqdus(c) − sldus(c) − 5FBu /10000

ftmelt = sldus(c) + 273.15 − 5FBu /10000

2-1
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. 

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.

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%

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.

The fuel melting correlation in FRAPCON-3.5 has been revised to account for gadolinia additions 

within the fuel. Based on available literature regarding gadolinia-doped urania (Toyoaki (1999), 

the melting point is estimated with the correlation presented in Equation 2.1-6. 

(2.1-6)

Equation 2.1-6 relates the melting temperature (Tmelt in K) to the concentration of gadolinia 

(XGd2O3 in wt%).

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. Due to the limited 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.

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%

2 3 Tmelt = 3115 − 4.8* X Gd O

2-3
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. 

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

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

Toyoaki, Y., M. Hideharu, and Y. Masahiro. 1999. “Melting and Measurements of Solidification 

Point of UO2-Gd2O3 Solid Solutions Under Solar Furnace,” Journal of High Temperature Study, 

25 (2): 71-79

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.5, FRAPTRAN1.5, 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.

[ ] exp( E / RT) 2RT

YK E

K T

T exp( / T) 1

K exp( / T) FCP 2 D

3 D

2 2 2

2

1 + + − θ −

θ θ =

[K exp( E / RT)] 2

Y

2

K T

exp( / T) 1

K FENTHL 3 D

2

1 2 + + − θ −

θ =

2-4
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.5.

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.5, FRAPTRAN-1.5, and MATPRO. However, FRAPCON-3.5 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 

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

be in good agreement with the experimental data. 

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.

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

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.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.

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

2-7
Leibowitz, L., L.W. Mishler, 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.5, 

FRAPTRAN-1.5, 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 adopted by PNNL and used in the FRAPCON3.5 and FRAPTRAN-1.5 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.5 (Wiesenack and Tverberg, 2000). 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 (Lanning et al., 2005) and continue to be used in 

FRAPCON-3.5.

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; 

2-8
however, it is only used to calculate the melt temperature. Unlike the models used in 

FRAPCON-3.5 and FRAPTRAN-1.5, 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 also be observed in the second term of Equation 2.3-1, which 

represents the electronic contribution. 

( )( ) ( )( )

[ ( )]{ [( ) ] } 3 2

th

V

5.2997x10 T exp 13358/ T 1 0.169 13358/ T 2

A BT'' 1 3e

C

1 6.5 0.00469T' 1 D

D k

+ − + +

















+ + 















+ − − =

 −



2-9
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)

K = Ko • FD• FP• FM• FR





 



  −



 





+

+ = − T

16361

exp T

4.715x10

0.0375 2.165x10 T

1 K 2

9

o 4

























+





 



 = +

T

B

0.0643

B

1.09

1 T arctan

B

0.0643

B

1.09 FD

3.265

3.265





















− + − 



 





− = +

) 100

T 1200 1 exp(

1

3 0.019B

0.019B FP 1

1 (s 1)p

1 p FM

+ −

− =

) 80

T 900 1 exp(

0.2 FR 1 − +

= −

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

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

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.5 and FRAPTRAN-1.5 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 temperaturedependent 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

= 1 396 exp( Q / T)

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.5 with a data correlation presented by 

Massih (Massih et al., 1992) to include gadolinia content (Lanning et al, 2005). The modified 

NFI with gadolinia dependency used in FRAPCON-3.5 is presented in Equation 2.3-10.

(2.3-10)

2 4

95 CT DT

A BT f(Bu) g(Bu) h(T)

1 K − +

+ + + =

exp( F/ T) T

E

A BT f(Bu) (1 0.9exp( 0.04Bu))g(Bu) h(T)

1 K95 2 + − + + + − − =

exp( F/ T) T

E

A a gad BT f(Bu) (1 0.9exp( 0.04Bu))g(Bu) h(T)

1 K

2 −

+

+ ⋅ + + + − − =





2-11
where 

a = constant = 1.1599

gad = weight fraction of gadolinia

As applied in FRACON-3.5 and FRAPTRAN-1.5, the above models are adjusted for asfabricated 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 by PNNL 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 NFI model at high 

temperature. This adjustment by PNNL was based on the Institute for Transuranium Elements 

data from 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 PNNL 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 by PNNL (Lanning et al., 2005) 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)

K 1.0789*K *[d /{1.0 0.5(1 d)}] d = 95 + −

)

T

D exp( T

C

A(x) a gad B(x)T f(Bu) (1 0.9exp( 0.4Bu))g(Bu) h(T)

1 K

2

mod

95

+ −

+ ⋅ + + + − − =





2-12
T = temperature, K

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

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

= 1 396 exp( Q / T)

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.5 accounts for gad additions, but the

modified NFI model in FRAPTRAN-1.5 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.5

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. 

0.0132exp(0.00188T)

0.1148 a gad 1.1599x 0.0040B 2.475x10 (1 0.00333B)

0.92 K95 4

+

+ ⋅ + + + − Φ = −





2-13
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 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.

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

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

UO2.

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

unirradiated UO2.

MATPRO (UO2)

1

2

3

4

5

6

7

8

1 2 3 4 5 6 7 8

Measured Thermal Conductivity (W/m*K)

Predicted Thermal 

Conductivity (W/m*K)

Predicted = Measured Christensen Godfrey

Bates Gibby Weilbacher

Goldsmith and Douglas Hobson Ronchi

Modified NFI (UO2)

1

2

3

4

5

6

7

8

1 2 3 4 5 6 7 8

Measured Thermal Conductivity (W/m*K)

Predicted Thermal 

Conductivity (W/m*K)

Predicted = Measured Christensen Godfrey

Bates Gibby Weilbacher

Goldsmith and Douglas Hobson Ronchi

2-15
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. The thermal conductivity values are 

measured ex-reactor from small portions of fuel irradiated at various temperatures and burnup 

levels. 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.

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

UO2.

1

2

3

4

5

1 2 3 4 5

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

2-16
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. 

Modified NFI (UO2)

1

2

3

4

5

1 2 3 4 5

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

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.

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. 

Verification of UO2, UO2-Gd2O3, and MOX thermal conductivity models are further 

demonstrated by the comparisons to measured in-reactor fuel temperatures as a function of 

burnup presented in the FRAPCON-3.5 assessment documentation (Geelhood and Luscher, 

2014).

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

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

MOX modeled with 7 wt% of plutonium.

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

unirradiated MOX modeled with 7 wt% of plutonium.

MATPRO (MOX)

1

2

3

4

5

1 2 3 4 5

Measured Thermal Conductivity (W/m*K)

Predicted Thermal 

Conductivity (W/m*K)

Predicted = Measured 3 wt% Pu 15 wt% Pu

Duriez-Modified NFI (MOX)

1

2

3

4

5

1 2 3 4 5

Measured Thermal Conductivity (W/m*K)

Predicted Thermal 

Conductivity (W/m*K)

Predicted = Measured 3 wt% Pu 15 wt% Pu

2-19
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.

• Only applicable to low burnup fuel (<5 GWd/MTU) because thermal conductivity 

degradation is not included in the model. 

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, 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 

2-20
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 

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

Wiesenack, W. and T. Tverberg. “Thermal Performance of High Burnup Fuel – In-Pile 

Temperature Data and Analysis,” in Proceedings of the ANS International Topical Meeting on 

Light Water Reactor Fuel Performance, pp. 730-737. April 2000, Park City, Utah. 2000.

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.

2-21
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.

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

4 P = eσT

2-22
(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 ± 8.9 percent, while the error 

reported in MATPRO is ± 6.8 percent. 

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.

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.

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

e 0.78557 1.5263x10 T −5 = +

2-23
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.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.5, FRAPTRAN-1.5, 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 (Luscher and Geelhood, 2011). 

These updated constants are included in FRAPCON-3.5/FRAPTRAN-1.5 (FRAP). Table 2.5-1 

presents the original MATPRO constants and the updated constants for UO2, as well as the 

constants for PuO2.

L / L K T K K exp( E / kT) ∆ 0 = 1 − 2 + 3 − D

2-24
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

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.5/FRAPTRAN-1.5 (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.

∆L / L0 = ∆L / L0 (T m ) + 0.043⋅ FACMOT

[ ( )] m m

5 ∆L / L0 = ∆L / L0 (T m ) + 0.043 + 3.6x10 T − T + ∆T −

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

FRAPCON-3/FRAPTRAN with updated constants.

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

FRAPCON-3/FRAPTRAN.

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

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

2-26
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. Due to the limited amount of data, the same 

standard deviation is applied to the PuO2 data.

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.

Luscher, W.G., and K.J. Geelhood. 2014. Material Property Correlations: Comparisons between 

FRAPCON-3.4, FRAPTRAN 1.4, and MATPRO, NUREG/CR-7024, PNNL-19417, Pacific 

Northwest National Laboratory, Richland, WA. 

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

2-27
[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.5 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.5 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.5 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.5 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).

For FTEMP < 1000K (2.6-1)

For FTEMP ≥ 1000K (2.6-2)

In FRAPCON-3.5, 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.5 models 

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

et al., 1974).

(0.0015)RSNTR

L

L

m

 = − 

 



 ∆

(0.00285)RSNTR

L

L

m

 = − 

 



 ∆

100*RSNTR /(3.0* FDENS) L

L

m

 = 

 



 ∆

2-28
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)

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 have 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 

( )

(TSINT 1453)

22.2 100 DENS

L

L

m −

− −  = 

 



 ∆

( )

(TSINT 1453)

66.6 100 DENS

L

L

m −

− −  = 

 



 ∆

[ ( )] [ ( )] ( ) 3 FBU B 35 FBU B

m

e 2.0e

L

L

L

L − + − +  + + 

 



 ∆ = ∆

2-29
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.5 and MATPRO and incorporates Equations 

2.6-4, 2.6-5, and 2.6-6. 

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

-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

2-30
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.

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 

2-31
creep-induced elongation (FCREEP) and densification due to pressured sintering (FHOTPS) and 

irradiation (FUDENS) to calculate the overall dimensional changes in fuel.

The gaseous swelling correlation in FSWELL was not used in previous versions of FRAPCON-3 

because it significantly overpredicts swelling. However, a new gaseous swelling model was

devised for FRAPCON-3.5 after ramp tests suggested gaseous swelling may contribute to 

permanent cladding deformation in high burnup rods. As for solid swelling, a modified version 

of the solid swelling correlation presented in the MATPRO FSWELL subroutine was adopted for 

FRAPCON-3. Although this solid swelling model is still used in FRAPCON-3.5, an additional 

model and recommendations for modeling solid swelling in gadolinia-doped fuels is also 

provided. The FSWELL subroutine is only applicable to time scales on the order of minutes to 

hours and, therefore, not used in FRAPTRAN-1.5.

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)

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 it is assumed that the gas that causes 

swelling is 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.

Although previous releases of FRAPCON-3 have considered solid swelling, the effects of 

gaseous swelling have been ignored based on available model-to-data comparisons. Comparison

between FRAPCON-3 predictions and NRC-sponsored test rod results revealed that the 

MATPRO swelling model that included gaseous fission product buildup overpredicted test results 

at burnup levels up to 10GWd/MTU (Lanning and Bradley, 1984; Lanning et al., 1997). 

Consequently, it was determined that swelling could be modeled better by eliminating the 

contribution from gaseous fission products and increasing the contributions of solid fission 

S

29 SS 2.5x10 B − =

( ) [ ( )] [ ]

S

56 11.73 0.0162 2800 T 8.0x10 B Sg 8.8x10 2800 T e e B −27 − − − − = −

2-32
products from 0.669 percent (MATPRO and FRAPCON-2) to 0.77 percent (FRAPCON-3) per 10 

GWd/MTU. 

However, a gasesous swelling model has been introduced in FRAPCON-3.5. The FRAPCON3.5 model is based on data from Mogensen (Mogensen 1985) and was developed after ramp test 

results suggested gaseous swelling may influence permanent cladding hoop strain in high burnup 

rods. The linear strain is given as a function of temperature over the ranges given in equations 

2.7-3 and 2.7-4. These models are phased in between 40 and 50 GWd/MTU by applying a factor 

that varies linearly between 0 and 1 at 40 and 50 GWd/MTU, respectively.

5 2 4.55 10 4.37 10 − − = × − ×

∆ T

l

l (960° < T < 1370°C) (2.7-3)

5 2 4.05 10 7.40 10 − − = − × + ×

∆ T

l

l (1370° < T < 1832°C) (2.7-4)

As of FRAPCON-3.4, the solid swelling 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. This is the same assumption for solid swelling that is used in 

FRAPCON-3.5. The correlations used to calculate fuel swelling under these conditions are given 

in Equations 2.7-5 and 2.7-6, respectively.

For burnup < 80 GWd/MTU

(2.7-5)

For burnup > 80 GWd/MTU

(2.7-6)

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 

The variable “bus” in Equations 2.7-5 and 2.7-6 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.5 versions of FSWELL is presented 

in Figure 2.7-1 along with experimental data from immersion density measurements. Based on 

*(2.315 10 * 2.987 10 ) −23 −24 soldsw = bus × + sigswell ×

bus fdens*2.974x10 *(bu bul) 10 = −

*(3.211 10 *5.974 10 ) −23 −24 soldsw = bus × + sigswell ×

2-33
this figure, it can be seen that the FRAPCON-3.5 model predicts less swelling than the MATPRO 

model, which tends to overpredict the experimentally determined data points.

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.5 version of 

FSWELL. Data comparisons between FRAPCON-3.5 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.

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)

2-34
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.

The solid swelling model in FRAPCON-3.5 was also revised to model gadolinia-doped fuels. 

For gadolinia-doped fules, the swelling rate was reduced from 0.62 to 0.50% per 10 GWd/MTU 

based on swelling data from several integrated fuel assemblies irradiated in Halden (Colombier 

2010). This lower swelling rate for gadolinia-doped fuels is extended to higher burnups rather 

than modified to reflect any assumed changes in swelling behavior. It is also recommended that 

“zero” densification be assumed for gadolinia-doped fuels to reproduce experimentally observed 

swelling behavior.

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, 

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

17 GWd/MTU62 GWd/MTU62 GWd/MTU65 GWd/MTU68 GWd/MTU70 GWd/MTU75 GWd/MTU78 GWd/MTU90 GWd/MTU90 GWd/MTU99 GWd/MTU99 GWd/MTU2-35
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.

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.

Mogensen, M., C.T. Walker, I.L.F. Ray, and M. Coquerelle. 1985. “Local Fission Gas Release 

and Swelling in Water Reactor Fuel during Slow Power Transients” Journal of Nuclear Materials 

vol. 131 pp. 162-171

Colombier, V., M.I. Horvath, T. Tverberg, H. Jenssen. 2010 “Comparative Study of the Gadolinia 

Experiments IFA-515.10, IFA-636, IFA-681.1, and IFA-676.1” HWR-924, OECD Halden 

Reactor Project, Halden, Norway

2-36
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, 

Optimized ZIRLO, and M5. Various heat treatments can be accommodated by specifying the 

coldworked condition of the alloy. Examples of coldworked conditions for the different alloys 

are provided in Table 3.0-1. 

Table 3.0-1 Example heat treatments and coldworked conditions for different alloys. 

Alloy Heat Treatement Cold Worked 

Condition

Zr-4 Cold Worked, Stress Relief Annealed (CWSRA) 50%

Zirlo Cold Worked, Stress Relief Annealed (CWSRA) 50%

Optimized Zirlo Partially Recrystalized Annealed < 50%

M5 Recrystallized Annealed 0%

Zr-2 Recrystallized Annealed (typical Zr-2 condition) 0%

In addition, several FRAPTRAN-1.5 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.5, FRAPTRAN-1.5, and MATPRO is applicable to Zircaloy-2, Zircaloy-4, 

ZIRLO, Optimized ZIRLO, and M5 alloys. However, FRAPTRAN-1.5 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.5/FRAPTRAN-1.5. 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.5/FRAPTRAN-1.5 codes assume a 

standard error of 10 J/kg*K below 1090K.

3-1
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.5/FRAPTRAN-1.5 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.5/FRAPTRAN-1.5 codes.

In addition to Zircaloy alloys, FRAPTRAN-1.5 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, 

Optimized 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

3-2
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.

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, Optimized 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.

3-3
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.5, and FRAPTRAN-1.5. This correlation is applied to 

Zircaloy-2, -4, ZIRLO, Optimized ZIRLO, and M5. In addition, FRAPTRAN 1.4 includes a 

correlation to calculate the cladding thermal conductivity of ZrNb-1.

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 

2 5 2 9 3 k 7.51 2.09x10 T 1.45x10 T 7.67x10 T − − − = + − +

σk = 1.01

k = 36

σk = 5

3-4
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.5 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.5, and 

FRAPTRAN-1.5 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 

diffusivity data that were not included in development of the model. Thermal diffusivity data, 

such as that provided by Peggs and Bunnel. et al. (1983), 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.

The relative error determined by comparison between the data presented in both Figures 3.2-1 

and 3.2-2, respectively, and the CTHCON model described in Equation 3.2-1 is approximately 

9%.

 

0.000461843 * T k = 15.0636 e

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

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

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)

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)

3-6
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, 

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 

3-7
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(2001) and the FRAPCON-3.5/FRAPTRAN-1.5 codes. The MATPRO(2001)

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.5 and FRAPTRAN-1.5 is based on measurements obtained by 

Kingery (1954) from fully dense and porous (87 percent TD) zirconium oxides. 

3.3.1 Model Development and Comparisons

Oxide thermal conductivity models in MATPRO and FRAPCON-3.5/FRAPTRAN-1.5 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.5/FRAPTRAN-1.5 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.5/FRAPTRAN-1.5 is derived from data collected by Kingery

(1954). The FRAPCON-3.5/FRAPTRAN-1.5 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 

3-8
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.5 and FRAPTRAN-1.5 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.5/FRAPTRAN-1.5 (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 

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. 

K T10x81.1835.0 4

MATPRO

− +=

*T9599.1K 43.6(*T10x41.2( ))10x946.1*T10x 4 7 10

FRAP

− − − −= − −

3-9 
Figure 3.3-1 Oxide thermal conductivity measurements as a function of temperature compared with the MATPRO a

3.5/FRAPTRAN-1.5 (FRAP) cladding oxide thermal conductivity models. 

3-10 
3.3.2 Applicability and Uncertainty

The subroutine ZOTCON in FRAPCON-3.5/FRAPTRAN-1.4 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.5 and FRAPTRAN1.5. 

Liquid = 4.1K

3-11 
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). 

d < 3.88x10-6 (3.4-1) 

d > 3.88x10-6 (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 +=ε

1 =ε − d0.50808642.0





 



 − ε=ε

300

T1500

12 exp

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

1.0

0.9

0.8

0.7

0.6

0.5

0.4

0.3

0.2

0.1

0.0

0 20 40 60 80 100 120 

Emissivity (unitless)

Surface Layer Thickness (μm)

Murphy and Havelock Juenke and Sjodahl

Burgoyne and Garlick Uniform Oxide Burgoyne and Garlick Nodular Oxid

Burgoyne and Garlick with Crud Model Prediction < 1500 K

Juenke and Sjodah at 1575 K Model Prediction > 1500K

3-13 
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.5/FRAPTRAN-1.5 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.5/FRAPTRAN-1.5 codes by linear interpolation. 

In addition, the FRAPTAN-1.5 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 ±=σ 1.0





 



 − ±=σ

300

1500T

2 1.0 exp

3-14 
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 09822.177763.2 cos −





 



 



 





π − =ε +

3

33 10x

161

1083T 09822.177763.2 cos −





 



 



 





π − =ε +

6 2

11 10x04.1T10x7.9 − − =ε −

6 3

33 10x4.4T10x76.9 − − =ε −

0067.0

3

1

3

2

p 11 33 +ε+ε=ε

3-15 
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.5/FRAPTRAN-1.5 (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.5/FRAPTRAN-1.5 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.5/FRAPTRAN-1.5.

axial T 5 6 2.5060 10 4.4410 10 − − ε = − × + ×

diametral T 5 6 2.3730 10 6.7210 10 − − ε = − × + ×

axial T 3 6 8.300 10 9.70 10 − − ε = − × + ×

diametral T 3 6 6.800 10 9.70 10 − − ε = − × + ×

3-16
Figure 3.5-1 Model-to-data comparison for FRAPCON-3.5/FRAPTRAN-1.5 (FRAP) and 

MATPRO axial thermal expansion. Data from cladding tubes.

Figure 3.5-2 Model-to-data comparison for FRAPCON-3.5/FRAPTRAN-1.5 (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.

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

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

3-17
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.5/FRAPTRAN-1.5 (FRAP) codes to determine the Young’s modulus and the shear 

modulus, respectively. 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-18
3.6.1 Model Development and Comparisons

CELMOD

The CELMOD subcode used in MATPRO uses different forms of the same equations as those

used in FRAPCON-3.5/FRAPTRAN-1.5 (FRAP). The differences in equations 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.5 is coded to include this 

error, this value is not used. This standard error is not included in FRAPTRAN-1.5.

1 2 3

11 7 Y = (1.088x10 − 5.475x10 T + K + K )/ K

Y 9.21x10 4.05x10 T 10 7 = −

K = (6.61x10 + 5.912x10 T)∆ 11 8

1

K 2.6x10 C10

2 = −





 



 Φ = + − 3 25 10

K 0.88 0.12exp

3-19
The CELMOD subcode included in the FRAP codes differs only in form 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.5/FRAPTRAN-1.5)

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.5 and FRAPTRAN-1.5 use the same correlations to determine Young’s 

modulus in the subcode CELMOD. Although FRAPCON-3.5 includes the coding to return the 

standard error of the model, neither FRAPCON-3.5 nor FRAPTRAN-1.5 uses these errors. 

FRAPTRAN-1.5 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 celmod = (1.088x10 − 5.475x10 *ctemp + c *deloxy + c *cwkf)/ c

c (1.16x10 ctemp*1.037x10 )*5.7015 11 8

1 = +

c2 = 1.0

10 c3 = −2.6x10

( ( )) ( ) 25 25 c2 = 0.88* 1− exp − fnck /1x10 + exp − fnck /1x10

( )

3

1 2

10 7

K

4.04x10 2.168x10 T K K G − + + =

G 3.49x10 1.66x10 T 10 7 = −

3-20
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.5 and FRAPTRAN-1.5 use the same correlations to determine shear modulus 

in the subcode CSHEAR. Although FRAPCON-3.5 includes the coding to return the standard 

error of the model, neither FRAPCON-3.5 nor FRAPTRAN-1.5 returns the error. FRAPTRAN1.5 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.

K = (7.07x10 + 2.315x10 T)∆ 11 8

1

1 3 2

10 7 cshear = (4.04x10 − 2.168x10 *ctemp + c *deloxy + c *cwkf)/ c

11 8 c1 = 7.07x10 − ctemp* 2.315x10

c2 = 1.0

10 c3 = −0.867x10

( ( )) ( ) 25 25 c2 = 0.88* 1− exp − fnck /1x10 + exp − fnck /1x10

3-21
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 abofe 80 GPa (e.g., a quarter of Bolmaro’s data is collected at < 250°F). It should also 

be noted that, with exception to several data points near the lower bound from Busby (1966), the 

majority of data is from material in a recrystallized or annealed condition (i.e., zero coldwork). 

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-22
Figure 3.6-1 Predicted vs. measured values for Young’s modulus for various data sets. 

0

80

40

120

0 40 80 

Predicted Young's Modulus, GPa

Measured Young's Modulus, GPa

Lower Bound Measured = Predicted Upper Bound

Armstrong and Brown (Zirconium) Padel and Groff (Zirconium) Busby (Zr-4)

Spasic (Zr-2) Mehan (Zr-2) Northwood (Zr-2)

Bolmaro and Povolo (Zr-4)

3-23 
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.5/FRAPTRAN-1.5. 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 Zircaloy-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-24
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.

A modified version of the CAGROW correlation is used in FRAPCON-3.5 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.5 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.5 because 

significant irradiation axial growth is not expected under the short time scales considered by 

FRAPTRAN-1.5.

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.5 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.5. 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)

L / L A[exp(240.8 / T)]( t) (1 f )(1 0.02CW) z

1/ 2 ∆ = Φ − +

21 0.845 ax = 2.18x10 *Φ −

3-25
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.

It should be noted that these growth correlations should not be applied above the fluence levels of 

the data from which they are derived as different growth behavior may occur (e.g., accelerated 

growth at high fluence). Furthermore, it should be noted that the growth models are only 

applicable to fuel rods and should not be applied to other components that may be under an 

applied stress (e.g., PWR guide tubes or BWR tie rods).

21 0.81787 ax = 7.013×10 Φ −

25 0.98239 ax = 9.7893×10 Φ −

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

∆L/L (%) is 0.11.

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

of ∆L/L (%) is 0.07.

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

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

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.

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.

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

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

3-28
3.7.2 Applicability and Uncertainty

The CAGROW model in FRAPCON-3.5 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. The relative uncertainties are 22.3, 18.6, and 44.8%, 

respectively. These three alloys are typically used in PWRs, whereas recrystallized (RXA) Zr-2 

is typically used in BWRs. The absolute and relative uncertainties in growth strains for RXA Zr2 are 0.07 and 20.3%, respectively. The relative uncertainty is used for each material except 

ZIRLO, which incorporates the absolute uncertainty.

It should be noted that these growth correlations should not be applied above the fluence levels of 

the data as accelerated growth has been observed in some zirconium alloys at fluences greater 

than seen at current US NRC burnup limits. Furthermore, these models should only be used for 

fuel rods and not for other components that may have significantly different stress states (e.g., 

PWR guide tubes or BWR tie rods).

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.

3-29
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.8 Creep Rate (CREPR)

The subcode CREPR is used by FRAPCON-3.5 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.5, 

because creep occurs over much longer time intervals than those considered in FRAPTRAN-1.5

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.5 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 from 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 principal stresses could be included and account for the 

difference in creep behavior during tensile and compressive creep. 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 to account 

for differences that were not explained by thermal creep. 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 FRAPCON3.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)





 



 − 



 



 σ

ε = RT

Q exp E

a sinh

T

E A

n

i eff

th 

C f(T) 1 C2

eff

C

irr o ε = ⋅ φ ⋅σ ⋅

3-30
φ = 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.

Table 3.8-1 Parameters for FRAPCON-3.5 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.5 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)

th irr th irr ε = ε + ε +   

( ( )) 2.05

th irr

0.109

th irr

s

p 0.0216 2 tanh 35500 −+ + ε = ⋅ ε − ⋅ ε

(1 exp( 52 t)) t th irr th irr

s

H p ε = ε − − ⋅ ε ⋅ + ε ⋅ + +  

( ) 2 th irr th irr 1/

1/ 2

th irr

s

p

H exp 52 t

2 t

52

+ +

+ − ⋅ ε ⋅ + ε ⋅

⋅ ε ⋅ ε

ε =    

( )

2

i

2

o

2

o i

2

o

2

2 i

o o

2

i i

r r r

r

r r P P P r P r

−

− − +

σ =

3-31
(3.8-8)

(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 but not Optimized 

ZIRLO. Optimized ZIRLO is partially recrystallized, which can change creep characteristics. 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 SRA 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

o i

2

o

2

2 i

o o

2

i i

t r r

r

r r P P P r P r

−

− − −

σ =

2

i

2

o

2

o o

2

i i

l

r r

P r P r

−

− σ =

(( ) ( ) ( ) ) 2

r l

2

t r

2

eff 5 l t σ = 0. σ − σ + σ − σ + σ − σ

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

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 

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

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

3-33
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.5, and FRAPTRAN-1.5

codes. However, FRAPTRAN-1.5 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-34
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 L = ad

{ { [ ( )]}} 1 2 5 8 MH exp 2.6034x10 T 2.6394x10 T 4.3502x10 T 2.5621x10 − − − = + − + +

3-35
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.

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)

3-36
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.5, and FRAPTRAN1.5 codes use similar correlations to determine the gas thermal conductivity. However, 

FRAPCON-3.5 and FRAPTRAN-1.5 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.5 and FRAPTRAN-1.5 (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-1
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.5/FRAPTRAN-1.5 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)

( )

( ) ( )

1.08x10 *(T 273.15) 4.51x10 *(T 273.15) ]

1.009 * P / T / T 273.15 [17.6x10 5.87x10 * T 273.15

k P / T * 2.8516x10 9.424x10 *T 6.005x10 *T

7 2 11 3

2 2 4.2 4 5

8 10 14 2

− − −

− + + − +

= − + − +

− −

− −

− − −





6 1.45 5 9 1.3 k 4.44x10 *T 9.45x10 *(2.1668x10 *P / T) − − − = +

( ) ∑

∑ 



























+ − δ Ψ

=

=

n

i

n

j 1

i ij ij j

i i

mix

x 1 x

k x k

( )( )

( ) 















+

− − Ψ = φ + 2

i j

i j i j

ij ij M M

M M M 0.142M

1 2.41

1/ 2

j

2 i 3/

2 1/ 4

j

i

1/ 2

j

i

ij

M

M 2 1

M

M

k

k 1

















+

















































+

φ =

4-2
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.

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

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)

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

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)

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

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)

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

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

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)

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)

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

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

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)

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

4-7
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.

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.

4-8
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.

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.

4-9

Appendix A

Subroutine Source Codes

Appendix A Subroutine Source Codes

A.1 Uranium Dioxide/Mixed Oxide Properties

A.1.1 Fuel Melting Temperature (PHYPRP)

FRAPCON-3.5: PHYPRP

*deck phyprp

 subroutine phyprp(gadoln)

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-4.8*gadoln*100.0

A-1
 fdelta = 1.0e-10

 go to 110

100 c1 = comp

 ftmelt = sldus(c1)+273.15-5.0*fbu/10000.0-4.8*gadoln*100.0

 fdelta = liqdus(c1)-sldus(c1)-5.0*fbu/10000.0-4.8*gadoln*100.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.5: 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-2
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-3
c

A.1.2 Fuel Specific Heat Capacity (FCP) and Fuel Enthalpy (FENTHL)

FRAPCON-3.5: 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-4
 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.5: 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-5
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.5: 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-6
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.5: 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-7
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.5: FTHCON

*deck fthcon

A-8
 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-9
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-10
 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.5: 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-11
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-12
A.1.4 Fuel Emissivity (FEMISS)

FRAPCON-3.5: 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.5: 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

A-13
 return

c

 end

A.1.5 Fuel Thermal Expansion (FTHEXP)

FRAPCON-3.5: 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

A-14
c c1u, c2u, c3u, edu parameters updated by W.G. Luscher 10-17-07

c

 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.5: 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

A-15
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

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-16
A.1.6 Fuel Densification (FUDENS)

FRAPCON-3.5: 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

A-17
 go to 150

140 al2 = x

150 continue

 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.5: FSWELL

*deck fswell

 subroutine fswell (fdens,comp,bu,bul,ftemp,soldsw,gadoln)

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.987e-24)

 if(bu.ge.6912000.0) soldsw = bus*(3.211e-23+sigswell*5.974e-24)

 if(gadoln.gt.0.0) soldsw=bus*(1.867e-23+sigswell*2.315e-24)

 return

 end

A.2 Cladding Material Properties

A.2.1 Cladding Specific Heat (CCP)

FRAPCON-3.5: 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

A-18
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 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.5: 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

A-19
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

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

A-20
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

 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.5: (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

A-21
FRAPTRAN-1.5: (CTHCON)

*deck cthcon

 subroutine cthcon (ctemp, time, flux, coldw, ccon, cdkdt, 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 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)

A-22
 else

 ccon = 36.d0

 cdkdt = 0.d0

 endif

 return

c

 end

A.2.3 Cladding Oxide Thermal Conductivity (ZOTCON)

FRAPCON-3.5: 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.5: 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-23
A.2.4 Cladding Oxide Emissivity (ZOEMIS)

FRAPCON-3.5: ZOEMIS

*deck zoemis

 subroutine zoemis (ctmax,zroxid,emissv)

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 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.5: 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

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

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

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.5: CTHEXP

*deck cthexp

 subroutine cthexp (ctemp,cathex,cdthex)

c

c cthexp calculates axial and diametral thermal expansion of

A-25
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

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.5: CTHEXP

*deck cthexp

 subroutine cthexp (ctemp, cathex, cdthex, icm)

 implicit real (a-h,o-z)

A-26
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

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

A-27
 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

 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.5: 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

A-28
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 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

A-29
 if (celmod.lt.1.0) celmod = 1.0

 ucelmd = 6.4e+09

 return

 end

FRAPTRAN1.5: CELMOD

*deck celmod

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

A-30
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

 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

A-31
c

 end

FRAPCON-3.5: CSHEAR

*deck cshear

 function cshear (ctemp,fnck,cwkf,deloxy)

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 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

A-32
130 continue

 if (cshear.lt.1.0) cshear = 1.0

 ucsher = 9.0e+09

 return

 end

FRAPTRAN-1.5: CSHEAR

*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

A-33
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 = (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.5: 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:2 - Zircaloy 2

c 4 - Zircaloy 4

c 5 - M5

c 6 - ZIRLO

c 7 - Optimized ZIRLO

A-34
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

 ax2=7.013e-21*f2**0.81787

 elseif(icm.eq.6.or.icm.eq.7) 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

c Zircaloy-4

 if(siggro.gt.0.0.and.icm.eq.4) 

 &cagrow = cagrow*(1.0+siggro*0.223)

 if(siggro.lt.0.0.and.icm.eq.4)

 &cagrow = cagrow/(1.0-siggro*0.223)

c M5

 if(siggro.gt.0.0.and.icm.eq.5) 

 &cagrow = cagrow*(1.0+siggro*0.186)

 if(siggro.lt.0.0.and.icm.eq.5)

 &cagrow = cagrow/(1.0-siggro*0.186)

c ZIRLO and Optimized ZIRLO

 if(icm.ge.6) 

 &cagrow = max(0.0,cagrow+siggro*0.0005)

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.5: 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

A-35
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)

 &/(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*

A-36
 &(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

 end

A.2.9 Cladding Meyer Hardness (CMHARD)

FRAPCON-3.5: 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.5: 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

A-37
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

 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.5

*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.

A-38
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

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.

A-39
 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.5

*deck gthcon

 function gthcon (gmix,gtemp,gpres,gpthk)

 implicit real (a-h, o-z)

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)

A-40
 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)

 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))

 end

A-41




UNITED STATES

NUCLEAR REGULATORY COMMISSION

WASHINGTON, DC 20555-0001

------------------

OFFICIAL BUSINESS
NUREG/CR-7024, Rev. 1 Material Property Correlations: Comparisons between

FRAPCON-3.5, FRAPTRAN-1.5, and MATPRO