NUREG/IA-0453 

International 

 Agreement Report

Benchmarking of a Generic CANDU 

Reactor with PARCS, MCNP and RFSP

 

 

 

 

Prepared by: 

B. Arsenault, O. Shaikh, T. Downar*, D. Jabaay*, A. Ward*, Y. Xu*, 

AMEC-NSS, 

4th Floor, 700 University Ave, Toronto, Ontario M5G1X6, Canada

*University of Michigan, Nuclear Eng. and Rad. Science

1934 Cooley, Ann Arbor, MI 48109-2104, USA

K. Tien, NRC Project Manager 

 

Division of Systems Analysis

Office of Nuclear Regulatory Research

U.S. Nuclear Regulatory Commission

Washington, DC 20555-0001

Manuscript Completed: June 2015

Date Published: September 2015

Prepared as part of 

The Agreement on Research Participation and Technical Exchange

Under the Thermal-Hydraulic Code Applications and Maintenance Program (CAMP)

Published by

U.S. Nuclear Regulatory Commission
AVAILABILITY OF REFERENCE MATERIALS

IN NRC PUBLICATIONS

NRC Reference Material

As of November 1999, you may electronically access 

NUREG-series publications and other NRC records at 

NRC’s Library at 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 Publishing Office

Mail Stop IDCC

Washington, DC 20402-0001

Internet: bookstore.gpo.gov

Telephone: (202) 512-1800

Fax: (202) 512-2104

2. The National Technical Information Service

5301 Shawnee Rd., Alexandria, VA 22312-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 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

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 contractor prepared 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 under an international cooperative agreement for the exchange of technical 

information. 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.
NUREG/IA-0453 

International 

 Agreement Report

Benchmarking of a Generic CANDU 

Reactor with PARCS, MCNP and RFSP

 

 

 

 

Prepared by: 

B. Arsenault, O. Shaikh, T. Downar*, D. Jabaay*, A. Ward*, Y. Xu*, 

AMEC-NSS, 

4th Floor, 700 University Ave, Toronto, Ontario M5G1X6, Canada

*University of Michigan, Nuclear Eng. and Rad. Science

1934 Cooley, Ann Arbor, MI 48109-2104, USA

K. Tien, NRC Project Manager 

 

Division of Systems Analysis

Office of Nuclear Regulatory Research

U.S. Nuclear Regulatory Commission

Washington, DC 20555-0001

Manuscript Completed: June 2015

Date Published: September 2015

Prepared as part of 

The Agreement on Research Participation and Technical Exchange

Under the Thermal-Hydraulic Code Applications and Maintenance Program (CAMP)

Published by

U.S. Nuclear Regulatory Commission

iii

ABSTRACT

The purpose of this project is to benchmark a generic CANDU reactor with the computer codes 

PARCS, MCNP and RFSP. RFSP is a diffusion code used by the Canadian nuclear industry for 

fuel management and for safety analyses. A number of changes were made to PARCS and 

GenPMAXS, the utility used to generate the fuel tables for PARCS. The changes were required 

to allow the modeling of specific aspects related to CANDU reactors.

More specifically, the benchmarking exercise involves an assessment of the Finite Difference 

Method and the hybrid NEM/ANM solvers available in PARCS against the results of Monte 

Carlo simulations performed with MCNP. The results of the simulations calculated with the 

Finite Difference Methods in PARCS and RFSP are compared to determine possible areas of 

improvement in RFSP.

This work is an in-kind contribution and has been funded by Amec Foster Wheeler.

v

TABLE OF CONTENTS

Page

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

LIST OF FIGURES ................................................................................................................... vii

LIST OF TABLES......................................................................................................................ix

EXECUTIVE SUMMARY ...........................................................................................................xi

ACKNOWLEDGMENTS.......................................................................................................... xiii

ABBREVIATIONS ....................................................................................................................xv

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

2. GENERIC CANDU MODEL .............................................................................................. 2-1

2.1 Reactor Geometry.................................................................................................... 2-1

2.2 End Shield Region ................................................................................................... 2-2

2.3 Extrapolation Distances and Extrapolation Factors .................................................. 2-2

2.4 Fuel and System Parameters................................................................................... 2-2

2.5 Burnup Profile .......................................................................................................... 2-2

2.6 Vertical Control Devices ........................................................................................... 2-3

2.6.1 Liquid-Zone Controllers................................................................................ 2-3

2.6.2 Shutoff Rods................................................................................................ 2-3

2.6.3 Mechanical Control Absorbers..................................................................... 2-3

2.6.4 Adjuster Rods .............................................................................................. 2-4

3. INCREMENTAL CROSS SECTIONS ............................................................................... 3-1

3.1 Methodology ............................................................................................................ 3-1

4. LATTICE AND REFLECTOR PROPERTIES.................................................................... 4-1

4.1 RFSP Fuel and Reflector Tables.............................................................................. 4-1

4.2 PARCS Fuel and Reflector Tables ........................................................................... 4-1

4.2.1 Fuel Tables.................................................................................................. 4-1

4.2.2 Reflector Tables........................................................................................... 4-2

5. REACTOR-CORE CONFIGURATIONS............................................................................ 5-1

5.1 Reference Case – Configuration Number 1.............................................................. 5-1

5.2 Guide Tubes Inserted – Configuration Number 2 ..................................................... 5-1

5.3 Guide Tubes Inserted and LZCs at 50% - Configuration Number 3.......................... 5-1

5.4 Guide Tubes Inserted and LZCs at 100% - Configuration Number 4........................ 5-2

5.5 Adjusters In – Configuration Number 5..................................................................... 5-2

5.6 Asymmetric Case – Configuration Number 6............................................................ 5-2

6. CODE CHANGES REQUIRED IN PARCS AND GENPMAXS.......................................... 6-1

6.1 PARCS .................................................................................................................... 6-1

6.1.1 Modeling the Shell and Sub-shell................................................................. 6-1

6.1.2 Extrapolation Distances ............................................................................... 6-2

6.1.3 Moveable Devices........................................................................................ 6-3

6.2 GenPMAXS.............................................................................................................. 6-4

6.2.1 Processing the Cross-sections from WIMS-IST ........................................... 6-4

6.2.2 Process the Cross Sections from DRAGON-IST.......................................... 6-5

6.2.3 Generation of the Fuel Tables for PARCS ................................................... 6-5

7. SIMPLIFICATIONS USED IN THE FULL-CORE MODELS .............................................. 7-1

7.1 Fuel Density ............................................................................................................. 7-1
vi

7.2 Incremental Cross Sections in the Reflector Regions............................................... 7-1

7.3 Extent of the Models in the Axial and Radial Directions............................................ 7-1

7.4 Energy Discretization ............................................................................................... 7-2

8. RESULTS ......................................................................................................................... 8-1

8.1 Extrapolation Distances and Extrapolation Factors .................................................. 8-1

8.2 Reactivity Worths ..................................................................................................... 8-1

8.3 Power Profiles.......................................................................................................... 8-2

8.3.1 Channel Power Profiles (Thermal Extrapolation Distances Only)................. 8-2

8.3.2 Channel Power Profiles (Thermal and Epithermal Extrapolation

Distances).................................................................................................... 8-3

8.3.3 Bundle Power Profiles.................................................................................. 8-4

9. DISCUSSION.................................................................................................................... 9-1

10. CONCLUSIONS.............................................................................................................. 10-1

11. REFERENCES................................................................................................................ 11-1

APPENDIX A TABLES...........................................................................................................A-1

APPENDIX B FIGURES .........................................................................................................B-1

APPENDIX C PARCS USER’S MANUAL UPDATE ASSOCIATED WITH THE

BOUNDARY CONDITIONS .............................................................................C-1

APPENDIX D CODING ADDED INTO THE ANM DRIVER IN PARCS TO

ACCOMMODATE THE NEW BOUNDARY CONDITIONS ..............................D-1

APPENDIX E RESULTS OF THE TEST CASES USED TO VERIFY THE CODING

IMPLEMENTED IN PARCS FOR THE BOUNDARY CONDITIONS ................E-1
vii

LIST OF FIGURES

Page

Figure 1 Representation of a Fuel Channel.........................................................................B-1

Figure 2 Code Structure Implemented in the WIMS-to-PMAXS Module..............................B-1

Figure 3 Differences in Channel Power Calculated Between PARCS (FDM, 1-group

Boundary Condition) and RFSP-IST for the Reference Case

(PARCS-RFSP-IST)X100%/RFSP-IST .................................................................B-2

Figure 4 Differences in Channel Power Calculated Between PARCS (FDM, 1-group

Boundary Condition) and RFSP-IST for the LZCs at 100%

(PARCS-RFSP-IST)X100%/RFSP-IST .................................................................B-3

Figure 5 Differences in Channel Power Calculated Between PARCS (FDM, 1-group

Boundary Condition) and RFSP-IST for the Asymmetry (PARCS-RFSP-IST)

X100%/RFSP-IST.................................................................................................B-4

Figure 6 Channel Power Distribution Calculated with PARCS

(FDM, 1-group Boundary Condition) for the Asymmetric Case..............................B-5

Figure 7 Differences in Channel Power Calculated Between PARCS (FDM, 1-group

Boundary Condition) and MCNP for the Reference Case (PARCS-MCNP)

X100%/MCNP.......................................................................................................B-6

Figure 8 Differences in Channel Power Calculated Between PARCS (FDM, 1-group

Boundary Condition) and MCNP for the LZCs at 100% (PARCS-MCNP)

X100%/MCNP.......................................................................................................B-7

Figure 9 Differences in Channel Power Calculated Between PARCS (FDM, 1-group

Boundary Condition) and MCNP for the Asymmetry (PARCS-MCNP)

X100%/MCNP.......................................................................................................B-8

Figure 10 Differences in Channel Power Calculated Between PARCS

(Hybrid NEM/ANM, 1-group Boundary Condition) and MCNP for the

Reference Case (PARCS-MCNP)X100%/MCNP..................................................B-9

Figure 11 Differences in Channel Power Calculated Between PARCS

(Hybrid NEM/ANM,1-group Boundary Condition) and MCNP for the LZCs

at 100% (PARCS-MCNP)X100%/MCNP.............................................................B-10

Figure 12 Differences in Channel Power Calculated Between PARCS

(Hybrid NEM/ANM,1-group Boundary Condition) and MCNP for the

Asymmetry (PARCS-MCNP)X100%/MCNP........................................................B-11

Figure 13 Differences in Channel Power Calculated Between PARCS (FDM, 2-group 

Boundary Condition) and RFSP-IST (FDM, 1-group Boundary Condition)

for the Reference Case (PARCS-RFSP-IST)X100%/RFSP-IST..........................B-12

Figure 14 Differences in Channel Power Calculated Between PARCS

(FDM, 2-group Boundary Conditions) and RFSP-IST

(FDM, 1-group Boundary Condition) for the LZCs at 100%

(PARCS-RFSP-IST)X100%/RFSP-IST ...............................................................B-13

Figure 15 Differences in Channel Power Calculated Between PARCS

(FDM, 2-group Boundary Conditions) and RFSP-IST

(FDM, 1-group Boundary Condition) for the Asymmetry

(PARCS-RFSP-IST)X100%/RFSP-IST ...............................................................B-14

Figure 16 Differences in Channel Power Calculated Between PARCS

(FDM, 2-group Boundary Conditions) and MCNP for the Reference Case

(PARCS-MCNP) X100%/MCNP..........................................................................B-15

Figure 17 Differences in Channel Power Calculated Between PARCS

(FDM, 2-group Boundary Conditions) and MCNP for the LZCs at 100% 

(PARCS-MCNP)X100%/MCNP...........................................................................B-16
viii

Figure 18 Differences in Channel Power Calculated Between PARCS 

(FDM, 2-group Boundary Conditions) and MCNP for the Asymmetry 

(PARCS-MCNP)X100%/MCNP...........................................................................B-17

Figure 19 Differences in Channel Power Calculated Between PARCS

(Hybrid NEM/ANM, 2-group Boundary Conditions) and MCNP for the

Reference Case (PARCS-MCNP)X100%/MCNP................................................B-18

Figure 20 Differences in Channel Power Calculated Between PARCS 

(Hybrid NEM/ANM, 2-group Boundary Conditions) and MCNP for the

LZCs at 100% (PARCS-MCNP)X100%/MCNP ...................................................B-19

Figure 21 Differences in Channel Power Calculated Between PARCS

(Hybrid NEM/ANM, 2-group Boundary Conditions) and MCNP for the

Asymmetry (PARCS-MCNP)X100%/MCNP........................................................B-20

Figure 22 Differences in Bundle Power Calculated Between PARCS (FDM, 1-group

Boundary Condition) and MCNP for the Reference Case (PARCS-MCNP)

X100%/MCNP – Plane Number 1 .......................................................................B-21

Figure 23 Differences in Bundle Power Calculated Between PARCS (FDM, 1-group

Boundary Condition) and MCNP for the LZCs at 100% (PARCS-MCNP)

X100%/MCNP – Plane Number 1 .......................................................................B-22

Figure 24 Differences in Bundle Power Calculated Between PARCS (FDM, 1-group

Boundary Condition) and MCNP for the LZCs at 100% (PARCS-MCNP)

X100%/MCNP – Plane Number 4 .......................................................................B-23

Figure 25 Differences in Bundle Power Calculated Between PARCS (FDM, 2-group

Boundary Conditions) and MCNP for the Reference Case

(PARCS-MCNP)X100%/MCNP – Plane Number 1 .............................................B-24

Figure 26 Differences in Bundle Power Calculated Between PARCS (FDM, 2-group

Boundary Conditions) and MCNP for the LZCs at 100% (PARCS-MCNP)

X100%/MCNP – Plane Number 1 .......................................................................B-25

Figure 27 Differences in Bundle Power Calculated Between PARCS (FDM, 2-group

Boundary Conditions) and MCNP for the LZCs at 100% (PARCS-MCNP)

X100%/MCNP – Plane Number 4 .......................................................................B-26
ix

LIST OF TABLES

Page

Table 1 Summary of Reactor Physical Characteristics and Nominal Conditions..........A-1

Table 2 Moderator and Coolant Characteristics.................................................................A-2

Table 3 Nuclide Number Densities of Moderator and Coolant............................................A-2

Table 4 Keywords Used in the WIMS-IST TAPE16 File and Used by the GenPMAXS ......A-2

Table 5 Energy Group Structure Used in WIMS-IST and DRAGON-IST for the

Lattice Calculations...............................................................................................A-3

Table 6 Extrapolation Distances Calculated with MCNP and Used in PARCS (cm)...........A-4

Table 7 Extrapolation Factors Used in RFSP-IST..............................................................A-4

Table 8 Reactivity Worths Calculated with MCNP..............................................................A-4

Table 9 Reactivity Worths Calculated with RFSP-IST with Thermal Boundary

Conditions.............................................................................................................A-4

Table 10 Reactivity Worths with PARCS (Finite-Difference Method) with Thermal

Boundary Conditions Only ....................................................................................A-4

Table 11 Reactivity Worths with PARCS (Nodal Method) with Thermal Boundary

Conditions Only ....................................................................................................A-5

Table 12 Reactivity Worths with PARCS (Finite-Difference Method) with Fast and

Thermal Boundary Conditions...............................................................................A-5

Table 13 Reactivity Worths with PARCS (Nodal Method) with Fast and Thermal

Boundary Conditions.............................................................................................A-5

xi

EXECUTIVE SUMMARY

Amec Foster Wheeler and the U.S. NRC have an agreement for the assessment of the Code 

Application and Maintenance Program (CAMP). The agreement provides Amec Foster Wheeler

the access to the thermal-hydraulics and neutronics codes, and requires in exchange a 

contribution using the CAMP toolset.

Within this project, the calculations obtained using the finite-difference solver and the hybrid 

NEM/ANM solvers in PARCS are benchmarked against MCNP. Identical calculations are also 

performed with the standard toolset used in the Canadian industry, RFSP, using the finitedifference method. The comparison of the results between PARCS and RFSP is used to 

identify areas of possible improvement in RFSP.

A number of changes were required in the CAMP toolset and are documented in the report. 

The changes were required to allow the modeling of the characteristics of a CANDU reactor in 

PARCS. Changes were made to the GenPMAXS code to make sure that the macroscopic cross 

sections used in PARCS and RFSP are the same and can also be generated at any specific 

burnup.

The document also provides a description of the generic CANDU model used in the analysis. 

The generic CANDU model was built in PARCS, RFSP and MCNP for the benchmarking 

analysis. A number of core configurations are also identified to test the robustness of the 

solvers in PARCS and RFSP.

The extrapolation distances of the energy groups used in the diffusion codes were calculated 

with MCNP for a one-to-one comparison for the benchmarking exercise. The results of the 

analysis showed:

• Better agreement is obtained between PARCS and MCNP when the PARCS finitedifference solver is used as opposed to when the hybrid NEM/ANM method is used,

• The results of the calculations using the finite-difference solver in PARCS and RFSP are 

very similar to each other when the same boundary conditions are used. The similarities 

are in terms of bundle power, channel power and reactivity worths,

• The Canadian toolset could be improved by the use of epithermal and thermal boundary 

conditions calculated with MCNP, instead of using the thermal boundary condition for the 

epithermal and thermal groups, and

• The functionalities introduced in GenPMAXS and PARCS for the modeling of a CANDU 

reactor have been successfully tested and can be used in future analyses.

xiii

ACKNOWLEDGMENTS

The authors would like to acknowledge the financial support from Amec Foster Wheeler and the 

contributions from David Luxat and Julian Lebenhaft. Appreciation is also given to Jim Donnelly 

whose insights and experience contributed to the understanding of the results of this work and 

the resolution of issues identified during the course of the work.

xv

ABBREVIATIONS

AAs Adjuster Absorbers

ANM Analytic Nodal Method

CANDU CANada Deuterium Uranium reactor

CAMP Code Application and Maintenance Program

CR Control Rod

FDM Finite-Difference Method

GC Generic CANDU

GenPMAXS Generation of the Purdue Macroscopic XS Set

LZCs Liquid-Zone Controllers

MCAs Mechanical Control Absorbers

MCNP Monte Carlo N-Particle

NEM Nodal Expansion Model

PARCS Purdue Advanced Reactor Core Simulator 

PM Poison Moderator

RELAP Reactor Excursion and Leak Analysis Program

RFSP-IST Reactor Fueling Simulation Program – Industry Standard Toolset

SORs Shutoff Rods

U.S. NRC United States Nuclear Regulatory Commission

WIMS-IST Winfrith Improved Multi-group Scheme – Industry Standard Toolset

1-1

1. INTRODUCTION

Amec Foster Wheeler entered an agreement to participate into a work effort within the Code 

Applications and Maintenance Program (CAMP). The agreement is a work in-kind that provides 

Amec Foster Wheeler with the thermal-hydraulic and neutronic codes developed within CAMP. 

The thermal-hydraulic codes are the Reactor Excursion and Leak Analysis Program (RELAP5) 

[1] and the TRAC/RELAP Advanced Computational Engine (TRACE) [2]. The neutronic codes 

are the Purdue Advanced Reactor Core Simulator (PARCS) [3] and the Generation of the 

Purdue Macroscopic XS Set (GenPMAXS) [4].

The focus of this report is to document the work that has been done in the area of core 

neutronics to allow the specifics of CANada Deuterium Uranium (CANDU) reactors to be 

modeled in the GenPMAXS/PARCS suite of computer codes. The document provides a 

detailed description of the simplified model of a CANDU reactor containing a representation of 

the major features (dimensions, fuel type, reactivity devices) used in the study. Many minor 

support components and details of the fuel burnup distribution have been omitted from the 

model as they were not considered essential for this work.

The logic that was implemented in GenPMAXS/PARCS to model the specifics of CANDU 

reactors are documented in this report. Also, the Finite-Difference Method (FDM) and the hybrid

Nodal Expansion Model/Analytic Nodal Model (NEM/ANM) method are assessed to determine 

the most appropriate method of solving the neutronic flux in the CANDU reactor. The results of 

the simulations calculated with PARCS are then compared to MCNP5 [5]. The results of the 

simulations calculated with the FDM method in PARCS are also compared to the results 

calculated with the FDM method in the Reactor Fueling Simulation Program (RFSP) [6] code 

commonly used in the Canadian industry to determine the consistency between the two solvers.

An investigation on the use of appropriate extrapolation distances is also performed using the 

extrapolation distances calculated with MCNP. A recommendation is made for the most 

appropriate boundary conditions to be used in the diffusion codes in order to reduce the errors 

associated to the bundle-power and channel-power distributions.

2-1

2. GENERIC CANDU MODEL

A generic CANDU reactor that contains 480 channels with 13 bundles each was assembled and 

used to benchmark the PARCS and the RFSP computer codes for a number of core 

configurations. In order to isolate the principal differences in the neutronics codes PARCS and 

RFSP and to provide consistent comparisons with MCNP, only fresh fuel was used in the 

simulations and the moderator was un-poisoned. It is expected that any potential neutronics 

differences between PARCS and RFSP which might occur because of fuel burnup 

heterogeneities would be less than the differences which would be revealed from the 

heterogeneities which might occur because of reactor control devices, which were explicitely 

modeled in the generic CANDU.

The model is simpler than a CANDU model in that the device locations have been set to fit lines 

representing whole fuel bundles or half bundles, or aligned with mesh lines. To reduce the 

required number of meshes, many structural materials have been omitted. The positions of the 

devices are not exactly at the same relative positions as in a typical CANDU reactor. However, 

the position deviations are of the order of a few centimetres and are less than the physical 

dimensions of the devices. The geometrical dimensions of the Liquid-Zone Controllers (LZCs) 

have also been modified from the ones that are used in a typical CANDU. Consistent generic 

CANDU models have been developed for use with each of PARCS and RFSP.

2.1 Reactor Geometry

For both PARCS and RFSP, device locations are specified by defining boundaries in a XYZ 

mesh where the properties representing the impact of fuel and devices are placed. The 

incremental cross sections of the materials are calculated with a transport code over a lattice in 

the XY directions (28.575 cm) and one bundle length in the Z direction (49.53 cm). The devices 

are mapped to their respective locations in the core and the associated incremental cross 

sections are applied to their respective locations in the core. This is in contrast to geometry 

specifications in the MCNP computer code documented in Reference [5] in which the actual 

shapes and positions can be represented with the exact geometry. Since the Generic CANDU 

(GC) is not an actual reactor, the geometry specifications presented herein have been selected 

to place devices at half or integral fuel bundle positions in order to simplify the PARCS and 

RFSP modeling. This approach is taken to avoid use of complicated meshes that would be 

required to match the as-built locations of an actual reactor. It should be noted that this 

description slightly displaces some devices from their as-built locations in an actual reactor. 

However, the displacements are usually small, generally within the physical dimensions of the 

devices.

The GC has 480 fuel channels with 13 fuel bundles per channel. One half bundle is located 

outside the calandria vessel at each end of each fuel channel. The reactor length is thus 

defined by 12 bundles within the axial length of the calandria vessel. Furthermore, the reflector 

is not of uniform thickness in that, at either end of the cylindrical calandria vessel, there is a 

section of reduced diameter (the notch).

Based on the standard CANDU pitch (XY-plane) of 28.575 cm, a bundle length of 49.53 cm, a 

mesh structure of half-a-bundle meshes was defined for the X, Y and Z directions (14.2875 cm 

for the X and Y directions and 24.765 cm for the Z direction). All devices have been specified 

to fit within the uniform mesh lines. This causes small displacements of device positions from 

the locations in an operating reactor. 

Table 1 lists a summary of reactor physical characteristics and nominal operating conditions.
2-2

2.2 End Shield Region

The end shield region of the GC has also been simplified from the CANDU design. It consists of 

two stainless steel tubesheets, each 5 cm thick, penetrated by the fuel channel extensions to 

the outer tubesheets where the fueling machine connects, plus the endshield material which 

consists of a mixture of stainless steel ball bearings (60 vol%) and H2O (40 vol%) filling the 

space not occupied by the fuel channel extensions. The endshield region was not modeled in 

PARCS and RFSP, and instead an extrapolation distance was used in the axial direction. The 

same approach was used in the radial direction in PARCS and RFSP where an extrapolation 

distance was used at the boundary of the reflector region. 

2.3 Extrapolation Distances and Extrapolation Factors

The extrapolation distances were calculated with MCNP [5] using a full-core configuration of the 

generic CANDU reactor without any in-core components. The thermal neutron fluxes were 

tallied and the radial and axial extrapolation distances were calculated and defined as the 

distance from the edge of the core to where the neutron flux vanishes to zero. The two 

extrapolation distances calculated with MCNP in the radial and axial directions in units of cm 

were used directly in PARCS.

RFSP uses an extrapolation factor defined by Equation 2.1

ff = 2 × Σtttt2 × dd (2.1)

where

ff is the extrapolation factor

Σtttt2 is the macroscopic transport cross section for the thermal group (cm-1

)

dd is the extrapolation distance (cm) of the thermal group calculated with MCNP 

Two extrapolation factors were calculated for RFSP: one for the axial direction and one for the 

radial direction. The transport cross section used to calculate the extrapolation factor in the 

axial direction was the thermal transport cross section of the lattice. The transport cross section 

used to calculate the extrapolation factor in the radial direction was the thermal transport cross 

section of the reflector.

2.4 Fuel and System Parameters

A fuel site or fuel channel is made of a fuel bundle with the associated bundle end plates 

immersed in heavy water coolant and enclosed in a pressure tube. The pressure tube is 

located inside a calandria tube with a gap filled with Helium gas between the pressure and the 

calandria tube. 

Figure 1 represents a fuel channel with the different components.

Table 2 lists the coolant and moderator characteristics. Table 3 lists the nuclide number 

densities of the coolant and moderator.

2.5 Burnup Profile

Fresh Natural Uranium (NU) bundles were used in the whole reactor core at a zero burnup 

value. There was no burnup profile used in the reactor core models to eliminate an additional 

complexity when evaluating the differences in power profile and reactivity worth calculated with 

PARCS and RFSP.
2-3

2.6 Vertical Control Devices

Many of the control devices enter the reactor vertically from the top of the reactor. The following 

sections list the details of the Liquid Zone Controllers (LZCs), Shutoff Rods (SORs), Mechanical 

Control Absorbers (MCAs) and Adjuster Absorbers (AAs).

2.6.1 Liquid-Zone Controllers

For the purpose of spatial control, the reactor is divided into 14 zones. There are six zone 

controller units traversing the core, with four of them divided into two compartments and two of 

them divided into three, for a total of 14 compartments. Spatial power control is obtained by 

means of a column of light water (H2O) in each compartment. Bulk control is achieved by 

varying the water level in all compartments by the same proportion. Spatial power control is 

achieved by differential adjustment of the water level in individual compartments. The water 

level in each zone controller can be varied from 0 to 100 percent full as defined below.

In the Generic CANDU model, the integral geometry of an LZC is not modelled. The geometry 

has been simplified to a cylindrical Zircaloy tube with an inner diameter of 11.44 cm and an 

outer diameter of 11.68 cm. There are two lengths of LZCs, 642.9 cm or 700.09 cm, as 

appropriate for each zone control unit. The units extend downward from the XZ-plane. This 

means that the zone control guide tubes have the same lengths. The compartments are filled 

with H2O or He-filled (voided). The top compartment of each Zone Control Unit extends to the 

upper surface of the calandria vessel. The He-filled space above the surface of the H2O also 

extends to the upper surface of the calandria. In PARCS, there were no devices located in the

reflector region of the model, compare to the RFSP model where empty tubes were located in 

the reflector region above the top zone control compartments.

2.6.2 Shutoff Rods

The shutdown mechanism for the shutdown system consists of 32 vertical shutoff rods. In 

normal operation, these rods are fully withdrawn and available. 

Each shutoff rod unit comprises a shutoff rod suspended by a cable from a winch mechanism 

and a guide tube enclosing its path in the core region. The guide tube extension lies between 

the calandria top and the reactivity mechanism deck where the shutoff rod is normally poised. 

The guide tube is heavily perforated through the core region and is filled with D2O when the 

shutoff rod is poised. There are other springs and locators at the bottom end of the guide tube 

that were not included in the model.

The shutoff rod assembly consists of a Cd tube sandwiched between two stainless steel tubes, 

572 cm long. This rod is located inside a Zr-2 guide tube that has a 24.5% volume perforation.

All shutoff rods start from the same elevation and drop to the same elevation. In the Generic 

CANDU model, the lower edge of the shutoff rods will be poised at the level of the highest point 

of the reflector (0 cm reference height or 79.125 cm above the upper edge of the calandria 

tubes in the top row) and fall according to the representative SOR drop characteristic.

2.6.3 Mechanical Control Absorbers

There are four mechanical control absorbers. They are identical in structure and composition to 

the shutoff rods. Thus, the physical description of the shutoff rods also applies to the MCAs. 

The purpose of the MCAs is to provide coarse control override capability when a change in 
2-4

system reactivity exceeds the rate or depth capability of the LZCs.

2.6.4 Adjuster Rods

The generic CANDU has 24 adjuster rods. Adjuster rods are of three lengths according to their 

location in the core. Each adjuster unit comprises an element suspended by a stainless steel 

cable connected to a winch above the core, as for the SORs. The in-core section maintains the 

location within the confines of a perforated zirconium guide tube.

The adjuster rod assembly consists of a central rod (shim rod) of titanium inside a stainless 

steel tube. These two items together form the adjuster. They are not moved independently of 

each other. This rod is located inside a Zr-2 guide tube.
3-1

3. INCREMENTAL CROSS SECTIONS

The 2-group macroscopic incremental cross-sections for the LZCs, the AAs, the MCAs and the 

SORs have to be calculated for the 37-element fuel bundle in hot conditions. 

3.1 Methodology

The methodology for the calculation of the incremental cross-sections for the different types of 

reactivity devices, involves the use of DRAGON-IST [7] lattice cell code1

. A 3-Dimensional 

super-cell model is used in the DRAGON-IST calculations, representing two lattice cells 

containing a fuel bundle in each cell, and an orthogonal reactivity device between them. 

DRAGON-IST solves the neutron transport problem for this super-cell using the collisionprobability technique. The incremental cross-sections are calculated for the case where the 

super-cell is split into two regions: inner (containing the reactivity device) and outer (super-cell 

region not included in the inner region).

The homogenized cross sections are calculated for three configurations listed below.

1. No device at all,

2. Guide tube between the fuel channels,

3. Absorber and guide tubes between the fuel channels

The incremental cross sections are calculated from the cross sections of the inner region as 

follows:

Guide tube = set 2 – set 1

Absorber = set 3 – set 2

DRAGON-IST cannot calculate explicitly the cluster-of-tubes geometry in 3D calculations. 

Hence, the complex structure of the fuel and the reactivity devices must be cylindered, i.e. 

approximated with a cylinder of material obtained as a mixture of the materials used in the fuel 

or in the device. This process requires a two-step homogenization process that was applied 

using the following 2D models:

• 2D Fuel Cell Model, and

• 2D Super-Cell Model

In the 2D fuel cell model, the 37-element fuel bundle is located at the centre of a lattice cell, 

surrounded by a pressure tube, annulus gas, calandria tube, and moderator. The flux 

distribution obtained from the fuel neutron transport calculations for this geometry is used for 

material homogenization of the cell components. The cell materials are homogenized into three 

mixtures: homogenized fuel (includes the fuel, cladding and coolant), tube material (calandria 

tube, pressure tube and the CO2 gas between them) and moderator. The 89 energy-group 

macroscopic cross sections calculated for the homogenized materials are then used in the 2D 

super-cell model. All the incremental cross sections were based on the ENDF/B-VI data library2

and condensed to a two-energy group structure with a cutoff at 0.625 eV.

 1 Small differences are found between the lattice properties calculated with WIMS-IST (lattice properties used in 

RFSP) and DRAGON-IST. The assumption is that the use of the lattice properties calculated with WIMS-ISTor 

DRAGON-IST has negligible impact on incremental cross sections of the in-core structural materials. 2 The ENDF/B-VI data library used in the analysis was the version generated in Canada for the CANDU reactors.

4-1

4. LATTICE AND REFLECTOR PROPERTIES

Both PARCS and RFSP require a set of cross sections to model the 37-element Natural 

Uranium fuel bundles of the generic CANDU model under normal operating conditions. 

This section of the document provides the details required to understand how the fuel and 

reflector tables were generated for both computer codes. It is to be noted that only fuel tables 

and reflector tables were used in the full-core models with the use of extrapolation distances at 

the radial and axial boundary locations. No special fuel tables were used to represent the endshield regions.

4.1 RFSP Fuel and Reflector Tables

WIMS-IST [8] models a 2-D slice of a 37-element fuel bundle under normal operating conditions 

that are listed in 

Table 1. The model follows the guidelines for standard CANDU modelling with WIMS-IST with 

the only exception that fuel element endcaps were not included. The WIMS-IST lattice 

calculations were performed with the ENDF/B-VI nuclear data library1

.

A number of commands to the standard WIMS-IST utilities were used to perform postprocessing of the WIMS-IST output for RFSP as well as a binary file “Tape 16” that contains 

information saved after each cycle by the “cellav” and “region” keywords in the WIMS-IST input.

Only the compositions for the fresh core were used in the current analysis. Perturbations were 

also performed at different boron concentrations in the range of 0 ppm to 5 ppm, but not used in

the current simulations.

The reflector properties were calculated for the fresh bundles only with a range of boron 

concentration varying from 0 ppm to 5 ppm boron. However, the simulations performed in the 

current analysis were done only at a 0 ppm boron concentration. The reflector properties were 

also calculated at the conditions listed in 

Table 1.

4.2 PARCS Fuel and Reflector Tables

The purpose of this section is to document the use of GenPMAXS [4] to format the fuel cross 

sections and device incremental cross sections for the PARCS model of the generic CANDU 

reactor.

GenPMAXS can organize the cross sections such that the main fuel cross section file will also 

contain the device incremental cross sections for the LZCs (full and empty) as well as for the 

AAs, SORs, MCAs and guide tubes.

4.2.1 Fuel Tables

PARCS and GenPMAXS currently utilize a maximum of 12 state variables and only two 

variables were used for the fuel tables: Control Rod (CR) and Poison Moderator (PM). The CRs 

were used to define the incremental cross sections of the different devices and structural 

components and the PMs were used to calculate the lattice properties at different poison 

concentrations ranging from 0 ppm to 5 ppm boron in steps of 0.5 ppm. However, the current 

 1 The ENDF/B-VI data library used in the analysis was the version generated in Canada for the CANDU reactors.
4-2

simulations were all performed at a concentration of 0 ppm boron.

A total of 17 types of cross sections, or branches, were calculated for the fuel properties. The 

first branch corresponded to the lattice properties with 0 ppm boron and the values were 

tabulated for 43 different burnup values. Branches 2 to 7 corresponded to the incremental 

properties of the devices and were calculated at zero burnup only. The branches were: the 

guide tubes for the AAs, the AAs, the guide tubes for the SORs, the SORs, the empty LZCs, 

and the filled LZCs. Branches 8 to 17 corresponded to the fuel tables from 0.5 ppm boron to 5 

ppm boron in step of 0.5 ppm. The lattice properties at each boron concentration were 

tabulated at 43 different burnup values. However, only the properties for the fresh fuel and at a 

0 ppm boron concentration were used in the current simulations.

4.2.2 Reflector Tables

The moderator purity was used as a perturbation for the reflector tables. However, only the 

properties at 0 ppm boron were used for the current analysis. A total of 11 types of cross 

sections, or branches, were calculated for the reflector properties. The 11 branches were used 

to model the poison in the moderator from a 0 ppm concentration to a 5 ppm concentration in 

steps of 0.5 ppm. The reflector properties were calculated for the fresh fuel only as equilibrium 

burnup profiles were not used in the current analysis.

No specific branches were used to define structural materials or in-core reactivity components 

for the reflector region.
5-1

5. REACTOR-CORE CONFIGURATIONS

A few considerations were selected to determine a number of reactor-core configurations that 

could be used to assess the accuracy of the solvers in PARCS and RFSP. 

• The first consideration was to use a single burnup value in the whole core to simplify the 

problem. A fresh core was selected because no burnup calculations are required in the 

lattice cells, which introduces a peak at the centre of the core and challenges the 

neutron-flux solver. The other advantage is to avoid any errors or inaccuracies when the 

cross sections are interpolated at the given burnup value. 

• Full-core calculations without reactivity devices or structural components were required 

to assess the symmetry of the solution in both, PARCS and RFSP.

• Full-core calculations without reactivity devices and with the guide tubes were required 

to assess the impact of the structural materials on the numerical solution.

• Calculations with LZCs at 50% and 100% were required to calculate the reactivity worth 

of the light-water absorbers which is an important parameter in CANDU reactors. 

• Full-core calculations were used to assess the impact of the adjusters. 

• A side-to-side asymmetric case was also required to determine the accuracy of the flux 

solvers. The asymmetry was obtained by withdrawing all the adjusters from one side of 

the core.

Based on the considerations listed above, a total of six cases were established and used to 

assess differences in RFSP and PARCS predictions of reactivity worth and power profiles. The 

cases were run with the finite-difference method in RFSP and PARCS, and also solved with a 

nodal method in PARCS.

5.1 Reference Case – Configuration Number 1

The reference case consisted of a full-core with fresh bundles. There were no xenon cross 

sections used in the model and all the reactivity devices were withdrawn from the core. No 

structural materials were inserted in the core. All the parameters at the lattice-cell level were 

defined by the lattice code and not altered by the reactor-core calculations (i.e. fuel temperature,

coolant densities etc...).

5.2 Guide Tubes Inserted – Configuration Number 2

In this configuration, the guide tubes for the SORs, MCAs and AAs were inserted in the core. 

The empty LZCs were added as control devices by changing their bank percentage insertions in 

the core. In RFSP the incremental cross sections for the empty LZCs also extend to the 

reflector region. Incremental cross sections of reactivity devices and guide tubes were not 

applied in the reflector region in the PARCS full-core models.

5.3 Guide Tubes Inserted and LZCs at 50% - Configuration Number 3

In this configuration, the guide tubes were inserted in the core and the LZCs were set at an 

average level of 50%. The main objective was to determine the differences between PARCS 

and RFSP in reactivity worth and power profile.
5-2

5.4 Guide Tubes Inserted and LZCs at 100% - Configuration Number 4

Same configuration as in 5.3 with the only difference that the LZCs were at 100%. This case 

was used to analyze the linearity of the reactivity worth of the LZCs as a function of fill level.

5.5 Adjusters In – Configuration Number 5

In this configuration, all the guide tubes were inserted in the model, the LZCs were set at 0% 

and all the adjusters were inserted. The objective was to assess the accuracy of the flux 

solvers when the flux is flattened by in-core devices.

5.6 Asymmetric Case – Configuration Number 6

In this configuration all the guide tubes were in the core and the LZCs were at 0%. The 

adjusters were inserted on one side of the core and withdrawn on the other side.
6-1

6. CODE CHANGES REQUIRED IN PARCS AND GENPMAXS

A number of code changes were required to PARCS and GenPMAXS to introduce the 

capabilities required for modelling a CANDU reactor. The following sections provide the details

6.1 PARCS

PARCS is a three-dimensional reactor-core code which solves the steady-state and timedependent, multi-group neutron diffusion and SP3 transport equations in Cartesian, cylindrical 

and hexagonal geometries. PARCS is coupled directly to the thermal-hydraulics system codes 

RELAP5 and TRACE which provide the temperature and flow field information to PARCS during 

the transient calculations via the few group cross sections.

The major calculation features in PARCS include the ability to perform eigenvalue calculations, 

transient (kinetics) calculations, Xenon transient calculations, decay heat calculations, pinpower calculations, and adjoint calculations for commercial Light-Water Reactors (LWRs).

6.1.1 Modeling the Shell and Sub-shell

RFSP has the capability to model CANDU reactor geometry, while PARCS was originally 

developed as a light water analysis code with uniform geometry. There are several unique 

geometric features of a CANDU reactor, including a notch region, making the model an 

extruded cylinder. In PARCS, this extrusion leaves an annular region at the ends of the core 

that need to be filled with a material. 

6.1.1.1 Input Design

The design for the input is an addition to the planar region and assembly type input cards. The 

planar regions are then assembled together via the planar region assign card, creating the full 

geometry, and the assembly type card maps one full axial assembly to the location in the radial 

configuration location. Both of these inputs allow a dummy region to appear anywhere in the 

core. A zero anywhere in this mapped input is treated as a dummy, or null region in the 

computational domain and therefore is not part of the solution. The rest of the nodes are used 

to set up the linear system.

6.1.1.2 Functionality

The new input allows for irregular and non-uniform geometry to be modeled in PARCS. A zero

cross section identifier can be placed anywhere in any planar region, which will call a modified 

Nodal Expansion Method (NEM) or Finite Expansion Method (FEM) driver1

. The code zeroes 

out values in the specified dummy regions specified by the zeros. The drivers perform checks 

on each node to see if it is a dummy region or not, and solves the proper one or two nodes 

problems accordingly. If no dummy regions are present in the model, PARCS defaults to the 

original NEM sweeping methodology, leaving previous cases unaffected by these coding 

changes.

 1 What the PARCS documentation refers to as “Finite Expansion Method” is what more commonly referred to in the 

literature as “Finite Difference Method” (FDM). The latter expression is used in this report.
6-2

6.1.1.3 Verification

A total of eight verification test cases were built. Four cases were of 2-Dimensional geometries 

of 3X3 nodes, and the other four cases were 3-dimensional geometries of 9X9 nodes with 11 

axial planes. For each set of four cases, two used FDM solvers, and the other two used NEM 

solvers. For the 2-dimensional 3X3 cases, the two cases to compare were (1) the input 

geometry with 3X3 nodes with no dummy regions present, and (2) the input geometry with 9X9 

nodes with all nodes being dummy regions except the 3X3 nodes in the centre. This setup 

made the models exactly the same for comparison purposes. For the 9X9 cases in 3-

dimensions, the purpose was to verify that a similar solution is achieved when using black 

absorber cross sections. The geometry was smaller and similar to the generic CANDU model. 

In these cases, a notch was created using dummy regions and black absorber cross sections.

The comparison between the 3X3 dummy and no dummy region cases agree very well. The 

eigenvalues for the dummy and no dummy cases for both FDM and NEM solvers agreed to 

within a fraction of a mk. The flux values for the epithermal and thermal groups agreed to within 

0.1% for the NEM solver, and match to the third decimal place for the FDM cases. It should be 

noted that the neutronic meshing was increased to 20 in the x direction and 20 in the y direction 

to ensure an accurate answer.

The 3-dimensional case comparison between the black absorber cross sections that were used 

in the calculations and the dummy region compared fairly well. The eigenvalue difference was 

about 1.5mk. The axially integrated power shape shows approximately 1% difference between 

the cases. Both cases used a zero incoming current albedo boundary condition.

6.1.2 Extrapolation Distances

To achieve more accurate comparisons between PARCS and RFSP, additional boundary 

conditions were added to PARCS. RFSP uses one boundary condition, or extrapolation 

distance, and PARCS has zero incoming current, zero flux, and reflective boundary conditions. 

Since the boundary conditions are fairly rigid for PARCS and not easily comparable to RFSP, 

albedo boundary conditions were implemented in PARCS. This boundary condition was then 

used to provide the structure for implementing an albedo extrapolation distance equivalent to 

the extrapolation distance used in RFSP.

6.1.2.1 Input Design

Initially in PARCS, there were three albedo cards, ALBEDO_R, ALBEDO_ZT, ALBEDO_ZB, 

which have an energy group dependent input of net current divided by the surface flux. The 

same style of input was used for the creation of the albedo boundary conditions in Cartesian 

coordinates with some modifications. Albedo cards need to be specified for each face in 

Cartesian geometry; north, south, east, west, top and bottom. In these cards, more than one 

set of albedos can be specified. In the last line of each radial card, a map is used to apply the 

specified boundary condition set at that axial position. In the case of the axial cards, the card 

BC_ZT_MAP (or ZB) is to be used to map each radial assembly location to the corresponding 

boundary condition set. 

The same type of input style is used for the extrapolation distance cards. EXTRAP_D_R 

requires an extrapolation distance to be specified for each energy group, along with an 

extrapolation radius to be used in calculating the nodal albedos. The same type of map used 

for the radial albedo cards can be used here. EXTRAP_D_Z requires the extrapolation distance 

for each energy group. Appendix C lists the new cards that were added for the treatment of the 

boundary conditions in PARCS.
6-3

6.1.2.2 Computation

A new derivation of the NEM solution was needed and coded into PARCS to accommodate the 

new boundary conditions. This coding was added in the one node solver of the ANM driver in 

PARCS. Appendix D provides the derivation of the equations that were implemented in 

PARCS. 

For the boundary condition associated with the extrapolation distance, the extrapolation radius 

and distance were used with the node positions to calculate the node-wise extrapolation 

distance for each face. The distance was then used with the diffusion coefficient to calculate 

the equivalent albedo boundary condition.

The new albedo cards provide the first step towards albedo-equivalent extrapolation distance 

boundary conditions. They are consistent with the boundary conditions that are already in 

PARCS. When using the albedo cards, the albedo boundary conditions that are equivalent to 

the standard boundary conditions of reflective, zero flux, and zero incoming current give the 

same results as the standard ones. The albedo cards also allow for a wider selection of 

boundary conditions to be used when running a case. The cards also allow for every axial 

mesh to have different radial boundary conditions and every radial node to have different axial 

boundary conditions.

6.1.2.3 Verification

The verification process consisted of a series of simple test cases in one, two and three 

dimensions using the hybrid NEM/ANM and FDM solvers and comparing the results using the 

albedo boundary conditions to the results using the standard boundary conditions. The next set 

of verification cases was a 2X2 assembly using the extrapolation distances and an equivalent 

albedo model, a 5X5 case with extrapolation distances which was used to verify that the correct 

albedos were being used in each node, and a 1-dimensional case compared to an analytic 

calculation. Appendix E lists the results for each case.

For the FDM cases, the comparison between the two methods show that the results were the 

same within the sixth decimal. For the NEM cases, the results are very close in the 1-

dimensional cases, but differences in the results were found between the two solutions with the 

higher leakage in the next 2- and 3-dimensional cases. For the other cases, the neutronics 

meshing remained at one, while for the NEM cases in 2- and 3-dimensions the meshing was 

increased to 12 and 20 meshes, respectively. 

The agreement between the albedo boundary conditions and the legacy boundary conditions in 

PARCS, which have been verified and validated themselves, verifies the new method.

6.1.3 Moveable Devices

CANDU reactors have control devices that are located in vertical positions in the reactor and 

perpendicular to the fuel channels. The task was to use a simple way to input a large number of 

control devices in PARCS and be simple to use as well as versatile enough to accommodate 

devices that do not fall perfectly within the nodal boundaries. The objective was to use the 

existing linear control rod bank in PARCS to move the vertical devices in the core.

6.1.3.1 Input Design

The new lincr_bank input card requires the control device tip coordinates for the fully inserted 

and withdrawn positions, a type specifying whether the rod is to be treated like a circle, square, 

or rectangle, and the dimensions of the type. By using coordinate input positions, it provides the 
6-4

functionality necessary for control devices positioned in the core regardless of the nodal mesh. 

The type and dimensions complement the coordinate positions by providing the functionality to 

calculate the nodal rod fractions. These values are the only information required to convert the 

linear control rod into a custom control rod input, allowing the user to input multiple devices with 

ease.

6.1.3.2 Computation

The new cards were implemented as an extra step in the input processing of a PARCS input 

file. The names of the new cards are linear_cr and lincr_bank, which are analgous to the 

custom_cr and ccr_bank cards.

The linear control bank card is similar to other control rod inputs in PARCS. The specific nodal 

locations are controlled by the linear control rod input. The rod fraction to perturb each node 

with control rod branch cross sections is also controlled by the input. The insertion of the 

control rods is controlled by the bank_pos card. The control rod branch cross sections to be 

used for each bank are controlled by the crb_def and crb_type cards.

6.1.3.3 Verification

A total of 4 test cases with a sample CANDU reactor were used to verify the coding in PARCS. 

The first test case was a single control rod inserted from the top of the reactor and located at the 

centre. The second test case consisted of a single rod also inserted from the top, but with a 

displacement in the radial and axial directions. The third test case consisted of a control rod 

inserted from the bottom, but with a displacement in the radial and axial directions. The fourth 

test case consisted of a configuration with the four control devices with the additional guide 

tubes inserted from the top to the bottom.

Each of these test cases was run using the linear control rod input and the custom control rod

input. The eigenvalues and the control rod positions were independently checked using the 

linear control rod input parameters to make sure that the incremental cross sections were 

applied at the correct locations. There was an exact match between the location of the rods and 

the eigenvalues calculated using the linear and the cutom control cards. The lincr_bank card is 

considered verified and provides a more straight forward way to move the control devices in a 

CANDU reactor.

6.2 GenPMAXS

The GenPMAXS code was developed to process the cross sections generated by a lattice 

physics code into the PMAXs format that can be read by PARCS. The burnup-dependent 

macroscopic cross sections are read by PARCS from the PMAXS file prepared by the code 

GenPMAXS and then used in PARCS to calculate the node-wise flux and power distribution, 

which is then used to calculate the region-wise burnup for time advancing the macroscopic 

cross sections.

6.2.1 Processing the Cross-sections from WIMS-IST

Changes were introduced in GenPMAXS to read the cross sections generated with the WIMSIST lattice code and located in the TAPE16 file. 
6-5

Table 4 lists the keywords used in the TAPE16 file and processed by the GenPMAXS utility 

code.

The module WIMS-to-PMAXS was developed to read the cross sections from the TAPE16 files 

and convert into the PMAXS format. Figure 2 shows the structure implemented in the WIMS-toPMAXS module.

The WIMS-to-PMAXS module consists of two major subroutines, WIMS_scan, and WIMS_read. 

As only two group data will be processed in the current version of GenPMAXS, and all burnup, 

branch, and history information are given in the GenPMAXS input file, the dimensions of the 

data structure are already determined. In the subroutine WIMS_scan, the data structure for 

describing the contents of the TAPE16 files is assigned according to the information from the 

GenPMAXS input file. Memory is allocated according to the dimensions of the files. In the 

subroutine WIMS_read, the cross sections data are read by first checking the keywords listed in 

Table 4. The cross section data are stored in the PMAXS data structure and prepared for 

generating partial derivatives and printing into the PMAXS file.

6.2.2 Process the Cross Sections from DRAGON-IST

This task was to modify GenPMAXS to process the cross sections for the CANDU benchmark 

model. This included the generation of the cross sections for the reflector. GenPMAXS was 

also modified to merge the cross sections from WIMS-IST and the incremental cross sections 

from DRAGON-IST into one PMAXS file.

6.2.3 Generation of the Fuel Tables for PARCS

The modified version of GenPMAXS was used to read the WIMS-IST cross sections and the 

DRAGON-IST incremental cross sections to produce a single PMAXS cross section file for 

PARCS. The input follows the standard GenPMAXS format used in conjunction with the new 

%INCREME card as described in Section 6.2.2. 

7-1

7. SIMPLIFICATIONS USED IN THE FULL-CORE MODELS

The computer models of the generic CANDU reactor that were built in PARCS, RFSP-IST and 

MCNP were all made consistent within each other during the course of the analysis. A few 

assumptions were used and justified below. 

7.1 Fuel Density

The modelling of each pin with the specific burnup in MCNP forces the simplification of the 

model where the burnups need to be binned to limit the number of universes in the model. In 

order to isolate the principal differences in the neutronics codes and to provide consistent 

comparisons with MCNP, fuel burnup was not considered.

The lattice cell calculations performed in WIMS-IST used a fuel density of 10.51 g/cm3 and an 

endcap card was used to take into account the end caps and the coolant on both ends of the 

bundle. The lattice properties in WIMS-IST are calculated over a bundle length and the use of 

the endcap card scales down the fuel density by taking into account the length of the fuel stack 

over the length of the bundle. The isotopic concentrations that were calculated in WIMS-IST in 

units of aa aaaa

bb bb×cc with the associated density of 10.5038 ggcc 3 was used in MCNP.

7.2 Incremental Cross Sections in the Reflector Regions

CANDU reactors are characterized by vertical reactivity devices inserted from the top. To each 

device is associated a guide tube located from the top reflector region to the bottom reflector 

region. To simplify the models, the incremental cross sections of the guide tubes located at the 

bottom reflector region were omitted. The incremental cross sections of the devices located in 

the top reflector region were included in RFSP-IST and MCNP and were omitted in PARCS. 

The incremental cross sections of the devices located in the reflector region are calculated 

using a standard supercell model in DRAGON-IST that comprises the device surrounded by two 

fuel channels. The use of the standard supercell model in DRAGON-IST generates small 

incremental fission cross sections that are not fully representative of the physics phenomena in 

the reflector region. In MCNP this does not affect the physics because each device is modelled 

with its material composition and incremental cross sections are not required. In RFSP-IST, the 

fission cross sections are set to 0 in the reflector region. PARCS does not accept fission cross 

sections in the reflector region and as a direct consequence, no devices were introduced in the 

reflector region. 

The incremental cross sections associated with the guide tubes located in the top reflector 

region are very small because they consist of a zircaloy material filled with heavy water. 

However, the incremental cross sections associated with the empty liquid-zone controllers is not 

negligible because the device consists of a zircaloy tube filled with void. The absence of heavywater in the top LZC compartments located in the reflector region is expected to be noticeable 

at the top of the reactor. 

7.3 Extent of the Models in the Axial and Radial Directions

The MCNP computer code has few limitations in terms of modelling capabilities and can 

represent complex geometries where the diffusion approximation does not hold. The model in 

the axial direction included the end shield region. The model also included the light water 

surrounding the reactor in the radial direction. 
7-2

In the PARCS and RFSP-IST core models, the extent of the model in the axial direction stopped 

at the calandria side tubesheet where an extrapolation distance was applied. In the radial 

direction, the models stopped at the surface of the heavy-water reflector where an extrapolation 

distance was applied. Two sets of extrapolation distances were used in the PARCS and in the 

RFSP-IST models; one set in the axial direction and one set in the radial direction.

The impact of using extrapolation distances in the PARCS and in the RFSP-IST models should 

have a small impact on the flux profile and power distribution as long as the extrapolation

distances are calculated using the equivalent MCNP full-core model. 

7.4 Energy Discretization

An advantage of the MCNP computer code is the use of a continuous energy group and 

eliminates the need for self-shielding calculations. WIMS-IST and DRAGON-IST both require 

discrete energy group calculations and apply the 89-energy group structure using the energy 

cutoffs listed in Table 5 for the lattice calculations. For the supercell calculations, the energy 

group are condensed into a 33-energy group structure to optimize the computer resources to 

solve the 3-Dimensional supercell configurations. Once the supercell problem is solved in the 

33-energy-group structure, the cross sections are condensed into a 2-energy group strcture to 

define the thermal and epithermal incremental cross sections used in the full-core diffusion 

calculations. The energy cutoff used to define the energy group structure is at 0.625eV. 
8-1

8. RESULTS

The first subsection provides the extrapolation distances that were calculated with MCNP with 

the derived extrapolation factors for use with the P card of the *DATA GEOMETRY module in 

RFSP-IST. The next sections provide the results of the simulations with MCNP, PARCS and 

RFSP-IST. The results include the reactivity worths of the different in-core devices and the 

power profiles using two sets of boundary conditions. In the first set, the methodology was 

consistent with the approach used in RFSP-IST where the boundary condition for the thermal 

group was also applied to the epithermal group. In the second set, the extrapolation distances 

calculated for the thermal and for the epithermal groups were applied in PARCS for the thermal 

and epithermal groups, respectively.

8.1 Extrapolation Distances and Extrapolation Factors

The extrapolation distances were calculated with MCNP as described in Section 2.3 and are 

listed in Table 6.

The numerical values listed in Table 6 representing the average extrapolation distances were 

directly used in PARCS because the definition of the extrapolation distances used in PARCS is 

the same as what was calculated in MCNP. The extrapolation distances of the thermal group 

listed in Table 6 were converted into extrapolation factors using Equation 2.1. Table 7 lists the 

extrapolation factors that were used in RFSP-IST.

The boundary conditions listed in Table 6 and Table 7 were used in PARCS and RFSP-IST, 

respectively. The next sections provide the reactivity worths and the power profiles calculated 

with the diffusion codes.

8.2 Reactivity Worths

The reactivity worths of the devices were calculated in units of mk using Equation 8.1 below.

ρρdddddd = kkeeee

iiii −kkeeeeee

rr

kkeeeeee

iiii ×kkeeeeee

rr × 1000 (8.1)

Where

ρρdddddd is the reactivity worth in mk of the device,

kk eeee

ii is the k-effective value of the full-core model with the device inserted, and

kk eeee is the k-effective value of the reference full-core model with no devices in.

Table 8 lists the reactivity worths that were calculated with MCNP and used as a reference to 

benchmark the results calculated with PARCS and RFSP-IST. Table 9 and Table 10 list the 

reactivity worths that were calculated with PARCS and RFSP-IST with the finite-difference 

method using the thermal boundary conditions only. The results show that the two codes 

predict very similar results. The reactivity worth of the LZCs is overpredicted by 8.5% and 8.1% 

in RFSP-IST and PARCS, respectively. The worth of the adjusters is underpredicted by 1.7% 

and 3.0% in RFSP-IST and PARCS, respectively.

The results listed in Table 11 indicate that the use of the nodal method increases the 

discrepancy of the reactivity worth associated with the LZCs. The worth of the LZCs is

overpredicted by 12.01% compared to MCNP. However, the use of the nodal method 
8-2

decreases the discrepancy of the reactivity worth associated with the adjusters compared to 

MCNP. The FDM method predicts a difference of 2.97% compared to MCNP and the nodal

method predicts a difference of 1.49%. 

Table 12 and Table 13 list the results calculated with PARCS when the two-group boundary 

conditions are used in the FDM and nodal methods, respectively. The reactivity worths listed in 

the two tables show that the use of the 2-group boundary conditions has a negligible impact on 

the reactivity worths when compared to reactivity worths calculated with the thermal group 

boundary condition only.

8.3 Power Profiles

The power profiles that were calculated with MCNP for the different configurations were used as 

a benchmark to determine the accuracy of the solutions obtained with PARCS and RFSP-IST. 

The finite-difference solver and the hybrid NEM/ANM solver where tested in PARCS using two 

sets of extrapolation distances. In the first set, the thermal extrapolation distances were applied 

to the fast and thermal groups. In the second set, the thermal extrapolation distance and the 

epithermal extrapolation distance were applied to the thermal and fast-neutron groups, 

respectively.

The channel power distributions calculated with PARCS were benchmarked against MCNP for 

all the configurations. For the bundle power distributions, the planes located at the extremity of 

the reactor in the axial direction were benchmarked against MCNP. Those planes were 

selected because of their sensitivity to the extrapolation distances. Plane number 4 was also 

selected because of the presence of the LZCs.

All the figures included in this Section represent a front view of the reactor where columns 1 to 

24 represent the fuel columns from the left-hand side to the right-hand side of the reactor. 

Rows A to Y represent the fuel rows where row A is located at the top of the reactor and row Y 

is located at the bottome of the reactor1

.

The next sections present the results of the benchmarking exercise and cover the similarities 

and differences that were identified in channel-power profiles and bundle-power profiles, 

repectively. 

8.3.1 Channel Power Profiles (Thermal Extrapolation Distances Only)

The thermal extrapolation distances that were calculated in the radial and axial directions for the 

thermal group only and listed in Table 6 where used in PARCS for the thermal and epithermal 

groups. The extrapolation factors listed in Table 7 for the thermal group were used in RFSP-IST

for the thermal and epithermal groups.

8.3.1.1 Benchmarking of PARCS Against RFSP-IST

Comparison of the results calculated with PARCS and RFSP-IST with the finite-difference solver 

with symmetric configurations are listed in Figure 3 and Figure 4. The results show that PARCS 

and RFSP-IST agree very well for the bare lattice model where no devices are inserted in the 

core (Figure 3). Figure 4 shows differences in power distribution located mainly at the top of the 

reactor. This is due to the absence of incremental cross sections in PARCS at the top of the 

reactor in the reflector region. This becomes evident for the largest perturbation (LZC guide 

tubes) with differences of the order of 2.5% at the top of the reactor.

 1 Row I is not used in the model. The generic CANDU reactor has 24 rows in total.
8-3

Figure 5 shows a larger power error on the left-end side of the reactor for the asymmetric core 

configuration. The left-end side of the reactor is characterized by a lower power compared to 

the right-end side of the reactor as shown in Figure 6.

8.3.1.2 Benchmarking PARCS Finite-Difference Method Against MCNP

Figure 7 to Figure 9 show the differences in power profiles that were calculated between 

PARCS and MCNP. The results show a slight asymmetry from side-to-side that was due to the 

convergence in MCNP. Figure 8 also show larger errors at the top of the reactor due to the 

absence of incremental cross sections in the top-reflector region in the PARCS model. 

Figure 9 shows that the FDM method overpredicts the power in the high-power zone located on 

the right-end side of the Figure. Figure 5 also shows that RFSP-IST predicts a higher power 

compared to PARCS in the same region and RFSP-IST would give larger differences compared 

to MCNP for the asymmetric case shown in Figure 9.

8.3.1.3 Benchmarking PARCS Hybrid NEM/ANM Method Against MCNP

Figure 10 to Figure 12 show the differences that were calculated between the hybrid NEM/ANM 

method in PARCS and the power distribution calculated with MCNP. The results show that the 

hybrid NEM/ANM method overestimates the power at the edge of the reactor and 

underrestimates the power at the centre compared to the FDM method. 

8.3.2 Channel Power Profiles (Thermal and Epithermal Extrapolation Distances)

The thermal extrapolation distances were calculated in the radial and axial directions for the 

thermal group and epithermal groups and listed in Table 6 where used in PARCS for the 

thermal and epithermal groups, respectively. The extrapolation factors listed in Table 7 for the 

thermal group were used in RFSP-IST for the thermal and epithermal groups.

8.3.2.1 Benchmarking of PARCS Against RFSP-IST

Figure 13 to Figure 15 show the differences that were calculated between PARCS with both 

boundary conditions (thermal and epithermal) using the FDM method and RFSP-IST using the 

FDM method with the thermal-energy-group boundary condition. The figures show that the 

differences in channel-power distributions were smaller compared to the differences in channelpower distributions using the 1-energy-group boundary condition (Figure 3). 

Figure 14 and Figure 15 show that the differences calculated between the PARCS FDM 2-

energy-group boundary conditions and RFSP-IST FDM 1-energy-group boundary conditions are 

similar to the differences calculated between the PARCS FDM 1-energy-group boundary 

conditions and RFSP-IST FDM 1-energy-group boundary condition (Figure 4 and Figure 5).

8.3.2.2 Benchmarking PARCS Finite-Difference Method Against MCNP

Figure 16 to Figure 18 show the differences calculated between the PARCS 2-energy-group 

boundary conditions using the FDM method and MCNP. The figures show that the differences 

in channel-power distributions were all similar to the differences that were calculated using the 

1-energy-group boundary condition (Figure 7, Figure 8 and Figure 9).
8-4

8.3.2.3 Benchmarking PARCS Hybrid NEM/ANM Method Against MCNP

Figure 19 to Figure 21 show the differences in channel-power distributions that were calculated 

between PARCS with the hybrid NEM/ANM solver and MCNP. The results show that the use of 

the hybrid NEM/ANM method increases the power error at the edge of the reactor in the lowpower regions compared to the FDM method (Figure 7, Figure 8 and Figure 9).

8.3.3 Bundle Power Profiles

The differences in bundle-power profiles were calculated between PARCS and MCNP. Only the 

calculations performed in PARCS with the FDM method were benchmarked against MCNP 

because Section 8.3.2.3 showed that the hybrid NEM/ANM method tends to overestimate the 

power for the bundles located at the periphery of the core. 

Figure 22 to Figure 24 show the differences in bundle powers that were calculated between 

PARCS FDM and MCNP when the thermal boundary conditions listed in Table 6 were applied 

to the fast and thermal groups in PARCS. Figure 22 shows for Plane number 1, located close to 

the end-shield region, differences close to -20% with an average value of -10.0% for Reference 

case. Figure 23 and Figure 24 show the differences in bundle powers for planes number 1 and 

plane number 4 with the LZCs at 100%. Figure 24 shows that near the centre of the core, the 

differences are small compared to MCNP and increase for the bundles located at the periphery 

of the core.

Figure 25 shows the differences in bundle power that were calculated between PARCS and 

MCNP at plane number 1 when the thermal and epithermal boundary conditions listed in Table 

6 were used for the thermal and epithermal groups. The results listed in Figure 22 and Figure 

25 show that the differences in bundle powers went down with the use of the 2-energy-group 

boundary conditions compared to the use of the thermal boundary condition only in PARCS. 

The average error in Figure 25 was -4.8% compared to -10.0% in Figure 22.

Figure 26 and Figure 27 show the differences in bundle power that were calculated between 

PARCS and MCNP at planes number 1 and number 4, respectively when the thermal and 

epithermal boundary conditions listed in Table 6 were used for the thermal and epithermal 

groups. Comparison of the results to Figure 23 and Figure 24, respectively show that the use of 

the 2-energy-group boundary conditions has a negligible impact in the regions located away 

from the end shields.
9-1

9. DISCUSSION

An improvement could be introduced in the full-core model in PARCS for the use of incremental 

cross sections of the devices located in the reflector region. This aspect is particularly important 

for the presence of the empty LZCs at the top of the reactor that reduces the amount of heavy 

water available for the moderation. 

The incremental cross sections associated with ννΣ could be set to zero in the reflector region, or 

a supercell model that extends into the reflector region could be used to calculate an 

incremental set compatible with the full-core model. The use of incremental cross sections at 

the top-reflector region in PARCS would increase the consistency between all the codes.

The extrapolation distances calculated in MCNP were based on the bare lattice model that has 

no in-core devices. It would be more consistent to calculate the extrapolation distances for each 

configuration.

10-1

10.CONCLUSIONS

Based on the results of the analysis, a number of observations can be drawn.

• The reactivity worth of the devices calculated with the finite-difference method with 

PARCS is approximately the same compared to the reactivity worth of the devices 

calculated with RFSP-IST. The differences are larger for the LZCs when the nodal 

method is used in PARCS. However, the use of the nodal method reduces the 

discrepancy of the reactivity worth associated to the adjusters compared to the FDM.

• The channel-power profile calculated with PARCS for the reference case when the finitedifference method is used compares very well with the channel-power profile calculated 

with RFSP-IST. 

• When devices or structural materials are located in the core, the maximum differences in 

channel-power profile calculated between PARCS and RFSP-IST is approximately 2.5%

at the top of the reactor close to the reflector region. This was due to the absence of 

incremental cross sections for the LZCs located at the top of the reactor in the reflector 

region in the PARCS model.

• The use of the nodal method in PARCS tends to increase the power of the channels 

located at the periphery of the core and slightly decrease the power of the channels 

located in the core.

• The power of the bundles located at the end of the channels and calculated with the 

FDM method in PARCS and RFSP-IST using the thermal extrapolation distances only is 

higher compared to the power of the bundles calculated with MCNP at the same 

locations. 

• The use of the thermal and epithermal boundary conditions with the FDM method in 

PARCS reduces the average error in bundle powers in plane number 1 located close to 

the end shield region.

• It is recognized that the nodal method predicts larger differences in power distribution at 

the periphery of the core compared to the finite-difference method when benchmarked 

against MCNP. Further investigation is required to understand the causes of the 

discrepancies and identify the most appropriate options associated with each method.

11-1

11.REFERENCES

[1] Information Systems Laboratories, Inc., “RELAP/MOD3.3 CODE MANUAL”, Volumes I 

to VII, March 2006.

[2] U.S. NRC, “TRACE V5.0 Theory Manual – Field Equations, Solution Methods, and 

Physical Models”, U.S. Nuclear Regulatory Commission, 2010.

[3] T. Downar, Y.Xu and V. Seker, “PARCS 3.0 U.S. NRC Core Neutronics Simulator –

User Manual”, September 2009.

[4] Y. Xu and T. Downar, “GenPMAXS-V5 – Code for Generating the PARCS Cross Section 

Interface File PMAXS”, University of Michigan, Report UM-NERS-10-0004, December. 

2009.

[5] X-5 Monte Carlo Team, “MCNP – A General Monte Carlo N-Particle Transport Code, 

Version 5”, LA-UR-03-1987, April 24, 2003.

[6] B. Rouben, “RFSP-IST, The Industry Standard Tool Computer Program for CANDU 

Reactor Core Design and Analysis”, Proceedings of the 13th Pacific Basin Nuclear 

Conference, Shenzhen, China, October 21-25, 2002.

[7] E G. Marleau, A. Hébert and R. Roy, “A User Guide for DRAGON 3.04T”, IGE-175, Rev. 

5T, August 2005.

[8] J.D. Irish and S.R. Douglas, “Validation of WIMS-IST”, Proc. of 23rd Annual Conference 

of Canadian Nuclear Society, Toronto, Canada, June 2-5, 2005.

A-1

APPENDIX A TABLES

Table 1 Summary of Reactor Physical Characteristics and Nominal Conditions

Calandria vessel

Inner radius (cm) 414.3375

Material Stainless Steel

Thickness (cm) 3.2

Core outer radius (cm) (variable - max) 363

Thickness of Reflector (cm) (variable - min) 51

Axial Length (cm) 594.36

Subshell inner radius (cm) 371.425

Axial extension of the sub shell – each end (cm) 49.5

End Shield Region

Tubesheets (both)

Thickness (cm) 5.0

Material Stainless Steel

Shield

Thickness (cm) 100.0

Material Stainless Steel/H2O

Ratio 60% Stainless Steel and 40% H2O

Fuel Channels

Number 480

Pitch (cm) 28.575

Length (cm) 

(12 in-core bundles+2 calandria-side tubesheets+2 endshield regions) 804.36 cm

Devices

Vertical

Shutoff rods (SDS1) 32

Adjuster rods 24

Zone controller units

Compartments/Zones

2-high units

3-high units

6

14

4

2

Mechanical control absorbers 4

Nominal Conditions

Moderator Purity (at%) 99.97

Coolant Purity (at%) 99.11

Moderator poison (ppm Gd) 0

Moderator temperature (K) 339

Coolant temperature (K) 563

Fuel temperature (K) 1011
A-2

Table 2 Moderator and Coolant Characteristics

Coolant

Purity (at%) 99.11

Temp (K) 563

Density (g/cc) 0.80768

Moderator

Purity (at%) 99.97

Temp (K) 339

Density (g/cc) 1.08709

Table 3 Nuclide Number Densities of Moderator and Coolant

Coolant

Isotope Number Density (at/barn.cm)

DD2O 4.8280E-02

H1H2O 4.3355E-04

O16 2.4299E-02

Moderator

Isotope Number Density (at/barn.cm)

DD2O 6.5490E-02

H1H2O 1.9652E-05

O16 3.2677E-02

Table 4 Keywords Used in the WIMS-IST TAPE16 File and Used by the GenPMAXS

Keyword Purpose

BEGIN Beginning of a cross section data block

ABSORPTION Absorption cross section, including fission

NU-FISSION Nu-fission cross section

TOTAL-X Transport cross section

FLUX Assembly average flux

SCATTER Scattering cross section

FISSPECT Fission spectrum

EDGEFLUX Assembly surface flux for ADF
A-3

Table 5 Energy Group Structure Used in WIMS-IST and DRAGON-IST for the Lattice 

Calculations

The Upper Energy Limit of Group 1 is 1.0E+07 eV

Group Lower Energy 

Limit eV Group Lower Energy Limit eV Group Lower Energy Limit Group Lower Energy Limit

1 7.7880E+06 24 9.1188E+03 47 4.0000E+00 70 3.0000E-01

2 6.0653E+06 25 5.5308E+03 48 3.3000E+00 71 2.8000E-01

3 4.7237E+06 26 3.3546E+03 49 2.6000E+00 72 2.5000E-01

4 3.6788E+06 27 2.0347E+03 50 2.1000E+00 73 2.2000E-01

5 2.8650E+06 28 1.2341E+03 51 1.5000E+00 74 1.8000E-01

6 2.2313E+06 29 7.4852E+02 52 1.3000E+00 75 1.4000E-01

7 1.7377E+06 30 4.5400E+02 53 1.1500E+00 76 1.0000E-01

8 1.3534E+06 31 2.7536E+02 54 1.1230E+00 77 8.0000E-02

9 1.0540E+06 32 1.6702E+02 55 1.0970E+00 78 6.7000E-02

10 8.2085E+05 33 1.3007E+02 56 1.0710E+00 79 5.8000E-02

11 6.3928E+05 34 1.0130E+02 57 1.0450E+00 80 5.0000E-02

12 4.9787E+05 35 7.8893E+01 58 1.0200E+00 81 4.2000E-02

13 3.8774E+05 36 6.1442E+01 59 9.9600E-01 82 3.5000E-02

14 3.0197E+05 37 4.7851E+01 60 9.7200E-01 83 3.0000E-02

15 2.3518E+05 38 3.7267E+01 61 9.5000E-01 84 2.5000E-02

16 1.8316E+05 39 2.9023E+01 62 9.1000E-01 85 2.0000E-02

17 1.4264E+05 40 2.2603E+01 63 8.5000E-01 86 1.5000E-02

18 1.1109E+05 41 1.7603E+01 64 7.8000E-01 87 1.0000E-02

19 8.6517E+04 42 1.3710E+01 65 6.2500E-01 88 5.0000E-03

20 6.7379E+04 43 1.0677E+01 66 5.000E-01 89 2.0000E-04

21 4.0868E+04 44 8.3153E+00 67 4.0000E-01

22 2.4788E+04 45 6.4760E+00 68 3.5000E-01

23 1.5034E+04 46 5.0435E+00 69 3.2000E-01
A-4

Table 6 Extrapolation Distances Calculated with MCNP and Used in PARCS (cm)

Epithermal Group Thermal Group

Radial Direction 4.13021 ± 0.25 3.55150 ± 0.07

Axial Direction 9.58116 ± 0.22 5.21171 ± 0.15

Table 7 Extrapolation Factors Used in RFSP-IST

Extrapolation Factor1

Radial Direction 2.85925

Axial Direction 3.94223

Table 8 Reactivity Worths Calculated with MCNP

Configuration Reactivity Worth (mk)

Guide Tubes -2.37 ± 0.022

LZCs at 50% -2.98 ± 0.03

LZCs at 100% -5.08 ± 0.02

Adjusters -22.19 ± 0.02

Core Asymmetry -7.54 ± 0.02

Table 9 Reactivity Worths Calculated with RFSP-IST with Thermal Boundary Conditions

Configuration Reactivity

Worth (mk)

(RFSP-ISTMCNP)/MCNP 

(%)

Guide Tubes -2.20 -7.3

LZCs at 50% -3.15 5.7

LZCs at 100% -5.52 8.5

Adjusters -21.81 -1.7

Core Asymmetry -7.36 -2.5

Table 10 Reactivity Worths with PARCS (Finite-Difference Method) with Thermal 

Boundary Conditions Only

Configuration Reactivity 

Worth (mk)

(PARCSMCNP)/MCNP 

(%)

Guide Tubes -2.19 -7.9

LZCs at 50% -3.15 5.5

LZCs at 100% -5.49 8.1

Adjusters -21.53 -3.0

Core Asymmetry -7.30 -3.3

 1 The extrapolation factors were calculated using Equation 2.1 with a thermal transport cross section of

 3.7821X10-1

cm-1 in the axial direction (fuel) and 4.0259X10-1 cm-1 in the radial direction (reflector). 2 The +/- values represent the uncertainties associated with the stochastic calculations in MCNP
A-5

Table 11 Reactivity Worths with PARCS (Nodal Method) with Thermal Boundary

Conditions Only

Configuration Reactivity 

Worth (mk)

(PARCSMCNP)/MCNP 

(%)

Guide Tubes -2.18 -8.1

LZCs at 50% -3.25 9.0

LZCs at 100% -5.69 12.0

Adjusters -21.86 -1.5

Core Asymmetry -7.36 -2.5

Table 12 Reactivity Worths with PARCS (Finite-Difference Method) with Fast and

Thermal Boundary Conditions

Configuration Reactivity 

Worth (mk)

(PARCSMCNP)/MCNP 

(%)

Guide Tubes In -2.18 -8.0

LZCs at 50% -3.15 5.8

LZCs at 100% -5.51 8.4

Adjusters In -21.43 -3.4

Core Asymmetry -7.28 -3.5

Table 13 Reactivity Worths with PARCS (Nodal Method) with Fast and Thermal

Boundary Conditions

Configuration Reactivity 

Worth (mk)

(PARCSMCNP)/MCNP 

(%)

Guide Tubes In -2.18 -8.1

LZCs at 50% -3.26 9.2

LZCs at 100% -5.71 12.3

Adjusters In -21.77 -1.9

Core Asymmetry -7.34 -2.7

B-1

APPENDIX B FIGURES

Figure 1 Representation of a Fuel Channel

Figure 2 Code Structure Implemented in the WIMS-to-PMAXS Module

WIMS_Scan

Assign the data structure to describe 

contents of cross section file

Build PMAXS branch structure from 

input branches

Memory allocation

WIMS_Read

Check the keywords

Read XS data
B-2

Figure 3 Differences in Channel Power Calculated Between PARCS (FDM, 1-group

Boundary Condition) and RFSP-IST for the Reference Case (PARCS-RFSPIST)X100%/RFSP-IST
B-3

Figure 4 Differences in Channel Power Calculated Between PARCS (FDM, 1-group

Boundary Condition) and RFSP-IST for the LZCs at 100% (PARCS-RFSPIST)X100%/RFSP-IST
B-4

Figure 5 Differences in Channel Power Calculated Between PARCS (FDM, 1-group

Boundary Condition) and RFSP-IST for the Asymmetry (PARCS-RFSPIST)X100%/RFSP-IST
B-5

Figure 6 Channel Power Distribution Calculated with PARCS (FDM, 1-group Boundary 

Condition) for the Asymmetric Case1

 1 Sum of all channel powers is equal to 1000.
B-6

Figure 7 Differences in Channel Power Calculated Between PARCS (FDM, 1-group

Boundary Condition) and MCNP for the Reference Case (PARCSMCNP)X100%/MCNP
B-7

Figure 8 Differences in Channel Power Calculated Between PARCS (FDM, 1-group

Boundary Condition) and MCNP for the LZCs at 100% (PARCSMCNP)X100%/MCNP
B-8

Figure 9 Differences in Channel Power Calculated Between PARCS (FDM, 1-group

Boundary Condition) and MCNP for the Asymmetry (PARCS-MCNP)

X100%/MCNP
B-9

Figure 10 Differences in Channel Power Calculated Between PARCS (Hybrid NEM/ANM, 

1-group Boundary Condition) and MCNP for the Reference Case (PARCSMCNP)X100%/MCNP
B-10

Figure 11 Differences in Channel Power Calculated Between PARCS (Hybrid NEM/ANM,

1-group Boundary Condition) and MCNP for the LZCs at 100% (PARCSMCNP)X100%/MCNP
B-11

Figure 12 Differences in Channel Power Calculated Between PARCS (Hybrid NEM/ANM,

1-group Boundary Condition) and MCNP for the Asymmetry (PARCSMCNP)X100%/MCNP
B-12

Figure 13 Differences in Channel Power Calculated Between PARCS (FDM, 2-group

Boundary Condition) and RFSP-IST (FDM, 1-group Boundary Condition) for

the Reference Case (PARCS-RFSP-IST)X100%/RFSP-IST
B-13

Figure 14 Differences in Channel Power Calculated Between PARCS (FDM, 2-group

Boundary Conditions) and RFSP-IST (FDM, 1-group Boundary Condition)

for the LZCs at 100% (PARCS-RFSP-IST)X100%/RFSP-IST
B-14

Figure 15 Differences in Channel Power Calculated Between PARCS (FDM, 2-group

Boundary Conditions) and RFSP-IST (FDM, 1-group Boundary Condition)

for the Asymmetry (PARCS-RFSP-IST)X100%/RFSP-IST
B-15

Figure 16 Differences in Channel Power Calculated Between PARCS (FDM, 2-group

Boundary Conditions) and MCNP for the Reference Case (PARCSMCNP)X100%/MCNP
B-16

Figure 17 Differences in Channel Power Calculated Between PARCS (FDM, 2-group

Boundary Conditions) and MCNP for the LZCs at 100% (PARCSMCNP)X100%/MCNP
B-17

Figure 18 Differences in Channel Power Calculated Between PARCS (FDM, 2-group

BoundaryConditions) and MCNP for the Asymmetry (PARCSMCNP)X100%/MCNP
B-18

Figure 19 Differences in Channel Power Calculated Between PARCS (Hybrid NEM/ANM, 

2-group Boundary Conditions) and MCNP for the Reference Case (PARCSMCNP)X100%/MCNP
B-19

Figure 20 Differences in Channel Power Calculated Between PARCS (Hybrid NEM/ANM, 

2-group Boundary Conditions) and MCNP for the LZCs at 100% (PARCSMCNP)X100%/MCNP
B-20

Figure 21 Differences in Channel Power Calculated Between PARCS (Hybrid NEM/ANM,

2-group Boundary Conditions) and MCNP for the Asymmetry (PARCSMCNP)X100%/MCNP
B-21

Figure 22 Differences in Bundle Power Calculated Between PARCS (FDM, 1-group

Boundary Condition) and MCNP for the Reference Case (PARCSMCNP)X100%/MCNP – Plane Number 1
B-22

Figure 23 Differences in Bundle Power Calculated Between PARCS (FDM, 1-group

Boundary Condition) and MCNP for the LZCs at 100% (PARCSMCNP)X100%/MCNP – Plane Number 1
B-23

Figure 24 Differences in Bundle Power Calculated Between PARCS (FDM, 1-group

Boundary Condition) and MCNP for the LZCs at 100% (PARCSMCNP)X100%/MCNP – Plane Number 4
B-24

Figure 25 Differences in Bundle Power Calculated Between PARCS (FDM, 2-

group Boundary Conditions) and MCNP for the Reference Case (PARCSMCNP)X100%/MCNP – Plane Number 1 
B-25

Figure 26 Differences in Bundle Power Calculated Between PARCS (FDM, 2-group

Boundary Conditions) and MCNP for the LZCs at 100% (PARCSMCNP)X100%/MCNP – Plane Number 1
B-26

Figure 27 Differences in Bundle Power Calculated Between PARCS (FDM, 2-group

Boundary Conditions) and MCNP for the LZCs at 100% (PARCSMCNP)X100%/MCNP – Plane Number 4
C-1

APPENDIX C PARCS USER’S MANUAL UPDATE ASSOCIATED

WITH THE BOUNDARY CONDITIONS

Card Type Field Default Description

ALBEDO_N Bcsetnum None The number of albedos to be input is specified 

by bcsetnum. The albedo is specified, for each 

bcsetnum, for each energy group for the (north, 

south, east, or west) face of the problem

domain. The last line of the input is a linear 

“map” of albedo set numbers specifying which 

albedo set will be used for that axial position.

In the example, three sets are input, one albedo 

for each of the two energy groups. In plane 1, 

set 1 is used, in plane 2 set 2, 3 set 3, 4 set 2, 

and 5 set 1.

albn(bcsetnum,ng)

irbcmap(nz)

ALBEDO_S bcsetnum None

albs(bcsetnum,ng)

irbcmap(nz)

ALBEDO_E Bcsetnum None

albe(bcsetnum,ng)

irbcmap(nz)

ALBEDO_W bcsetnum None

albw(bcsetnum,ng)

irbcmap(nz)

ALBEDO_* 3

0.0 0.0

0.33 0.4

0.5 0.5

1 2 3 2 1

ALBEDO_T bcsetnum None The number of albedos to be input is specified 

by bcsetnum. The albedo is specified, for each 

bcsetnum, for each energy group for the (top or 

bottom) face of the problem domain. 

In the example, two sets are input, one albedo 

for each of the two energy groups. The top and 

bottom albedos are to be used in conjunction 

with the BC_Z*_MAP card(s).

albt(bcsetnum,ng)

ALBEDO_B bcsetnum None

albb(bcsetnum,ng)

ALBEDO_T 2

0.0 0.0

100000.0 1000000.0

BC_ZT_MAP izbcmap(nfa) nfa*1 The mapped input should be the same 

dimensions as the planar region, or the number 

of radial assemblies. The albedo set number 

input at each assembly position will determine 

which set of albedos to use.

In the example, node 1 uses set 1, node 2 uses 

set 2, node 3 uses set 1, nodes 4-6 use 2, node 7 

uses 1, node 8 uses 2, and node 9 uses 1. 

BC_ZB_MAP izbcmap(nfa) nfa*1

BC_ZT_MAP

1 2 1

2 2 2

1 2 1
C-2

EXTRAP_D_R bcsetnum None The number of extrapolation distances to be 

input is specified by bcsetnum. The extrap. 

dist. is specified, for each bcsetnum, for each 

energy group. The last line of the input is a 

linear “map” of the set numbers specifying 

which set will be used for that axial position.

The example has 2 sets of extrapolation 

distances, specifying distances for each energy 

group and a geometric radius to use for each 

set. Plane 1 uses set 1, 2-3 use 2, and 4 uses 1.

extrapdistr

(bcsetnum,ng)

extraprad(bcsetnum)

Irbcmap(nz)

EXTRAP_D_R 2

2.0 2.0 28.0

2.5 2.7 29.0

1 2 2 1

EXTRAP_D_Z bcsetnum None The number of extrapolation distances to be 

input is specified by bcsetnum. The set is 

specified, for each bcsetnum, for each energy 

group of the problem domain. Should be used 

in conjunction with the BC_Z*_MAP card(s).

In the example, two sets are input, one extrap. 

dist. for each of the two energy groups.

extrapdistz

(bcsetnum,ng)

EXTRAP_D_Z 2

3.5 3.7

2.6 2.43
D-1

APPENDIX D1 CODING ADDED INTO THE ANM DRIVER IN PARCS

TO ACCOMMODATE THE NEW BOUNDARY

CONDITIONS

For NEM, we start with the transverse integrated 1D diffusion equation,

2

2 ' ' ,' '

' 1 ' 1

( ) ( ) () () () () ()

G G

g g tf

tg g tg g g g g g g g

g g

d x

D x x x Sx S x Lx

dx k

φ χ φ νφ φ = =

− +Σ − Σ − Σ = = − ∑ ∑

(1)

Define the vectors and matrices,

1 1

2 2 , ,

G G

J

J

J

J

φ

φ φ

φ

  

  

= =   

  

     

 

 

1

2 ,

G

S

S

S

S

 

 

=  

 

   





A T = 1 T S F

k − − χ

 

,

1

2

G

χ

χ

χ

χ

 

 

=  

 

   





1

2

f

f

fG

F

ν

ν

ν

  Σ

  Σ   =  

    Σ  





, 

1,1 1,2 1,

2,1 2,2 2,

,1 ,2 ,

,

G

G

G G G G

S

  ΣΣ Σ   ΣΣ Σ =

ΣΣ Σ  





 



T

1

2

0 0

0 0 ,

0 0

t

t

tG

  Σ

  Σ =

  Σ





 

 and 

D

1

2

0 0

0 0 ,

0 0 G

D

D

D

 

 

=

 





 



Then, diffusion equation (1) can be rewritten in matrix form,

− D

2

2 () () () d x A x Sx

dx φ φ + =   

(2)

Define non-dimensional variable , 2 2

rl rl x xx xx h

h h ξ + − =− =

The Diffusion equation becomes

D −

2

2 2 () () () d A S

h d φξ φξ ξ ξ

+ =   

(3)

Approximate the transverse leakage and transient fix source with 2nd order polynomial

( )

2

,

0

( ) g gi i

i

S sP ξ ξ == ∑

( ) ( ) ( ) ( ) ( ) 2 3 4 2

01 2 3 4

11 1 1, , (3 1), (5 3 ), (35 30 3) 22 8

PP P P P ξ ξξ ξ ξ ξ ξ ξ ξ ξ ξ = = = − = − = −+

 1 Equations listed in the Appendix that constitute the coding implemented in PARCS was provided by University of 

Michigan
D-2

( ) ( ) ( ) ( ) ( ) ( ) ( ) 0 1 2 3 24 3 P P P P PP P ' 0, ' 1, ' 3 , ' 5 1, ' 7 3 ξ ξ ξξ ξ ξ ξ ξξ = = = =+ =+

( ) ( ) ( ) ( ) ( ) ( ) ( ) 0 1 2 03 4 2 P P P PP P P " 0, " 0, " 3 , " 15 , " 35 10 ξ ξ ξ ξ ξξξ ξ = = = = = +

0

c s s =

1

1 ( )( 2 )( ) ( )( 2 )( ) c r cc lc l l cc rc r s Gh s s h h h h s s h h h h − = − + +−− + +    

( )

2 1

2 ( )( ) ( )( ) c r cc l l cc r s G h s sh h s sh h − = − ++− +    

2( )( )( ) c lc rc l r G h hh hh h h = + + ++

The Diffusion equation can be rewritten as,

D − ( ) 2 2

2 2

0

() () i i

i

d A sP

h d φξ φξ ξ ξ =

+ = ∑   

(4)

Approximate the solution with 4th order Legendre polynomials

( )

4

0

( ) i i

i

φξ ξ a P

=

= ∑  

(5)

The 5G coefficients of Legendre polynomials will be determined by 5G constraints. 

For the boundary node, these constraints are average flux, 0th (nodal neutron balance), 1st and 

2nd order WRE and the boundary condition. 

a0 = φ

 

Define: 

β D = h

     ,

Κ D =

1

2 h

−

      are diagonal matrices.

2nd order WRE: 4

1

35

a = Κ  ( Aa s 2 − 2 )  

0th order WRE: 2 4 3 10 a a + =Κ   ( Aa s 0 − 0 )  

a I 2 = +Κ 21 2  ( ) 1A −

Κ (7 72 Aa s s 0 − +0 2 )   

1st order WRE: 3 − +Κ 15a

 Aa1 = Κ  s1



For right boundary,

r r J = α  r

df r φ



Where, 

r α , r

df are diagonal matrices for albedo and discontinuity factors respectively. r γ r = α r

df

 

r D J = −  ( ) 4

1

i i ' 1 i

a P

h β =

∑ = −  (aa a a 12 3 4 +++ 3 6 10 )   
D-3

( )

4 4

01234

0 0

1 r

i i i

i i

φ aP a a a a a a

= =

= = = ++++ ∑ ∑    

− β ( 12 3 4 3 6 10 ) r aa a a +++ = α    r

df (aaa aa 01234 ++++ )  

3 a C = −  0 0 a C −  1 1 a C −  2 2 a C −  4 4 a



C 0

r γ =

6 β r + γ , C 1

β =

r + γ

6 β r + γ , C 2

3 β =

r + γ

6 β r + γ , C 4

10 β =

r + γ

6 β r + γ

a C 1 = − 15  ( 1 + Κ ) 1

A C 15 −

0a C 0 +  2a C 2 +  ( 4a4 ) − Κ  s1    



where CC0, CC1, CC2 aaaa CC4 are diagonal matrices

All the coefficients are then known to calculate the current. The current is then used to update 

the diffusion coefficient. All of this was coded into a separate option for the one-node NEM 

solver in PARCS.

E-1

APPENDIX E RESULTS OF THE TEST CASES USED TO VERIFY

THE CODING IMPLEMENTED IN PARCS FOR THE

BOUNDARY CONDITIONS

Boundary Cond. Albedo PARCS BC Albedo PARCS BC Albedo PARCS BC Albedo PARCS BC Albedo PARCS BC Albedo PARCS BC

Reflective 1.121065 1.121065 1.121065 1.121065 1.121065 1.121065 1.121065 1.121065 1.121065 1.121065 1.121065 1.121065

Zero Incoming J 0.304938 0.304938 0.386610 0.386610 0.186000 0.186026 0.204721 0.204721 0.091930 0.091946 0.130416 0.130416

Zero Flux 0.359818 0.359818 0.304938 0.304938 0.149503 0.149509 0.165962 0.165962 0.117134 0.117188 0.103049 0.103049

NEM

1-D 2-D 3-D

k-eff

FDM NEM FDM NEM FDM

K-eff

Group 1 6.622E+16 6.622E+16 6.622E+16 6.622E+16

Flux 6.622E+16 6.622E+16 6.622E+16 6.622E+16

Group 2 1.991E+16 1.991E+16 1.991E+16 1.991E+16

Flux 1.991E+16 1.991E+16 1.991E+16 1.991E+16

Two by Two Case

Albedo Extrapolation Distance

0.044034 0.044034

Distance Radius

Group 1 3.0 10.0

Group 2 3.0 10.0

Extrapolation Distance

0.681223

0.4418445

Two by Two Case Boundary Conditions

Albedo

x albedos G1 1 2 3 4 5 x albedos G2 1 2 3 4 5

1 1.090996 - - - 1.090996 1 0.707585 - - - 0.707585

2 0.211285 - - - 0.211285 2 0.137033 - - - 0.137033

3 0.170306 - - - 0.170306 3 0.110455 - - - 0.110455

4 0.211285 - - - 0.211285 4 0.137033 - - - 0.137033

5 1.090996 - - - 1.090996 5 0.707585 - - - 0.707585

y albedos G1 1 2 3 4 5 y albedos G2 1 2 3 4 5

1 1.090996 0.211285 0.170306 0.211285 1.090996 1 0.707585 0.137033 0.110455 0.137033 0.707585

2 - - - - - 2 - - - - -

3 - - - - - 3 - - - - -

4 - - - - - 4 - - - - -

5 1.090996 0.211285 0.170306 0.211285 1.090996 5 0.707585 0.137033 0.110455 0.137033 0.707585

Five by Five Case Calculated From Input

x albedos G1 1 2 3 4 5 x albedos G2 1 2 3 4 5

1 1.090996 - - - 1.090996 1 0.707585 - - - 0.707585

2 0.211285 - - - 0.211285 2 0.137033 - - - 0.137033

3 0.170306 - - - 0.170306 3 0.110455 - - - 0.110455

4 0.211285 - - - 0.211285 4 0.137033 - - - 0.137033

5 1.090996 - - - 1.090996 5 0.707585 - - - 0.707585

y albedos G1 1 2 3 4 5 y albedos G2 1 2 3 4 5

1 1.090996 0.211285 0.170306 0.211285 1.090996 1 0.707585 0.137033 0.110455 0.137033 0.707585

2 - - - - - 2 - - - - -

3 - - - - - 3 - - - - -

4 - - - - - 4 - - - - -

5 1.090996 0.211285 0.170306 0.211285 1.090996 5 0.707585 0.137033 0.110455 0.137033 0.707585

Five by Five Case From PARCS Debug Run





NUREG/IA-0453 Benchmarking of a Generic CANDU Reactor with PARCS, MCNP and RFSP