Reactor Core Methods

Yousry Azmy

•

Enrico Sartori

Nuclear Computational Science A Century in Review

13

Prof. Yousry Azmy North Carolina State University Department of Nuclear Engineering Raleigh, NC 27695 1110 Burlington Engineering Labs USA [email protected]

Enrico Sartori Organisation for Economic Co-operation and Development (OECD) 12 bd. des Iles 92130 Issy-les-Moulineaux France [email protected]

ISBN 978-90-481-3410-6 e-ISBN 978-90-481-3411-3 DOI 10.1007/978-90-481-3411-3 Springer Dordrecht Heidelberg London New York Library of Congress Control Number: 2009944067 Mathematics Subject Classification (2010): 82D75, 65C05 c Springer Science+Business Media B.V. 2010 No part of this work may be reproduced, stored in a retrieval system, or transmitted in any form or by any means, electronic, mechanical, photocopying, microfilming, recording or otherwise, without written permission from the Publisher, with the exception of any material supplied specifically for the purpose of being entered and executed on a computer system, for exclusive use by the purchaser of the work. Cover design: deblik Printed on acid-free paper Springer is part of Springer Science+Business Media (www.springer.com)

Ely Gelbard November 6, 1924–April 18, 2002

Preface

Scheduled on the heels of the atomic century, the American Nuclear Society’s international topical meeting on Mathematics and Computation seemed like an opportune moment in time to capture accomplishments in this area during the first half-century of nuclear engineering. Held in a semi-secluded part of the city of Gatlinburg, Tennessee, April 6–10, 2003, this gathering of prominent experts in the field and young professionals embarking on exciting careers in what promises to develop into a nuclear renaissance turned out to be the perfect venue for such a review. The conference was co-sponsored by three divisions of the American Nuclear Society, namely the Mathematics and Computation Division, the Reactor Physics Division, and the Radiation Protection and Shielding Division. The Technical Program of the conference revolved around the theme of its title, Nuclear Mathematical & Computational Sciences: A Century in Review, A Century Anew. The Anew component comprised contributed papers organized in 25 regular and special sessions on a broad variety of topics, plus a poster session and a panel session. The Review component of the conference comprised the lecture series that grew into this book. As Technical Program Chair (YYA) and Assistant General Chair (ES) of the conference, we decided to break with the traditional format of plenary sessions standard in technical meetings and organize a lecture series that takes stock of the state of the art in nuclear computational science at the turn of a new century. Thus the concept of the lecture series that led to the chapters of this book was born. One of the first experts we solicited to present a lecture in the series was the late Dr. Ely Gelbard of Argonne National Laboratory at the time. In his gentle, but firm and persuasive manner, he declined preferring instead to participate as co-organizer of the lecture series. We jumped on the opportunity recognizing his long-standing, distinguished, and generous contributions to many subareas in nuclear computational science, and his many years of service in the field positioned him well to know the major areas to cover in the lectures and to nominate world-renowned lecturers. In short order the three of us came up with a slate of topics and a corresponding list of lecturers. The response of the nominated lecturers was supportive and enthusiastic, and by mid Fall 2001 what has later become known as the Gelbard Lecture Series was fully conceived, and a tentative idea of ultimately documenting the lecture contents in book chapters was initiated. Our charge to the invited lecturers was to provide an overview of the assigned topic aiming primarily at breadth of coverage,

vii

viii

Preface

with a sharp focus on its mathematical and computational aspects. Specifically we requested that each author provide a historical perspective of the conception of their topic as a major area of research in nuclear computational science, and to identify landmarks for the evolution of the topic through the end of the twentieth century. We further requested that the lecturers delineate the current state of the art in their assigned topic and to project into the future by exposing perceived challenges and opportunities for advancing the frontier of knowledge. Our renowned lecturers did not disappoint and the lecture series was a smashing success, thanks to their dedicated effort and professionalism. The lectures, scheduled to open each half-day of the conference, were well attended, with conference participants packing the lecture hall on a consistent basis. Perhaps the only sour note that tainted the lecture series was the passing on April 18, 2002, of Dr. Ely Gelbard whose contributions to the success of the lecture series, and ultimately to the publication of this book, cannot be overstated. This great loss to the field of nuclear computational science overshadowed the conference leading to various observances of this sad event. The conference banquet included a memorial celebrating Dr. Gelbard’s life and his significant contributions to nuclear computational science, and the lecture series was named after him in recognition of his involvement that propelled the series to success. Later, the contributing authors to this book agreed to dedicate it to the memory of Dr. Ely Gelbard. Unfortunately death struck again with the passing of Dr. Richard Hwang on December 20, 2007, shortly after he completed the final revisions to his chapter appearing in this book. We are grateful for Richard’s contribution to the success of the lecture series, for the chapter he composed in this book, and for his dedication to his research over the past 5 decades. While the original list of topics envisioned in our early planning of the lecture series has not changed, the reader will notice a few differences between the lectures lineup and the chapters herein. First, Dr. Dan Cacuci who, for unforeseen circumstances, was unable to deliver his lecture on Sensitivity and Uncertainty Analysis at the conference has graciously composed the corresponding chapter for this book. Second, Dr. Kord Smith who presented the lecture on Reactor Physics at the conference apologized from composing the corresponding book chapter due to increased job-related responsibilities. We are grateful to Dr. Robert Roy for accepting to undertake such burden and for the excellent job he did in composing his chapter on Reactor Core Methods. Lastly, in composing Chapter 7, Elliott Whitesides recruited Mike Westfall and Calvin Hopper to help with the composition. This book would not have been possible without the support and active involvement of many people over the span of 6 years. Most of all we wish to thank the authors who willingly and cheerfully accepted this additional burden to their normally hectic schedules. We are confident that the benefit to the field of nuclear computational science and the gratitude of its practitioners, especially the young scientists who will carry the torch into the future, will reward the authors’ perseverance and patience during this long an arduous journey. We are grateful to Argonne National Laboratory’s Dr. Roger Blomquist for composing the memorials to Ely Gelbard and Richard Hwang, and for reviewing the final version of Richard’s

Preface

ix

Chapter 5. The support and encouragement of Bernadette Kirk, Director of Oak Ridge National Laboratory’s Radiation Safety Information Computational Center (RSICC) and General Chair of the Gatlinburg conference, was invaluable to the completion of this project. The technical help by Alice Rice of RSICC with bringing together the pieces of this book into a single volume is greatly appreciated. In addition, we wish to acknowledge the tacit approval and support of our respective institutions, The Pennsylvania State University and North Carolina State University (YYA), and the Nuclear Energy Agency of the Organisation for Economic Cooperation and Development (ES). June 2009

Yousry Y. Azmy Enrico Sartori

Obituary Composed by Dr. Roger Blomquist for Dr. Ely Meyer Gelbard

Ely Gelbard was born in New York City on November 6, 1924. He was the son of immigrants. His undergraduate work was at the City Colleges of New York and after World War II he earned his Ph.D. in physics from the University of Chicago. During the war, he served in the US Army Air Corps as a radar technician. He was a Senior Scientist at Argonne National Laboratory and a Fellow of the American Nuclear Society. Ely started his postgraduate career when the use of digital computers to solve the neutron balance equations for fission reactor core design and analysis was just starting to receive wide application. At Bettis (1954–1972), he participated in the efforts that put the numerical methods for the solution of the finite difference form of the neutron transport equation on a firm mathematical basis, and he devised several approximation schemes that were suitable for numerical methods and also developed efficient algorithms for their solution. While at Bettis, he earned international stature in the field, authoring important papers in many variants of the solution procedures (spherical harmonics, Sn , synthetic methods, and Monte Carlo), including the book, Monte Carlo Principles and Neutron Transport Problems, with J. Spanier. He was the first physicist at Bettis to attain the rank of Consulting Scientist, and earned the Atomic Energy Commission’s prestigious E. O. Lawrence Award. Since 1972, when Dr. Gelbard joined Argonne National Laboratory, fast reactors have been the focus of ANL’s reactor program, with its emphasis on more accurate computation of the neutron spectrum. His work in this area produced fundamental advances in the analysis of neutron streaming, collision probabilities, improvements in Monte Carlo methods, and neutron diffusion and transport within the nodal approximation. He also brought improved iterative solution strategies to bear on the equations of single-phase computational thermal-hydraulics analysis of passively safe metal-cooled reactor systems. He was consulted by many at ANL, at other labs, and at universities on a wide variety of technical issues, and invariably provided important insights. Ely’s sustained record of high productivity of the highest-quality technical work attracted a series of bright and vigorous visiting scholars and students whose participation magnified his work. He excelled at distilling complex technical issues to their essence, then performing the relevant mathematical analysis and, finally, computationally confirming the analysis. He was always careful, honest, and thoroughly

xi

xii

Obituary Composed by Dr. Roger Blomquist for Dr. Ely Meyer Gelbard

scrupulous in his work. He earned the ANS Special Award for Computer Methods for the Solution of Problems in Reactor Technology, the ANS Mathematics and Computations Division Distinguished Service Award, the ANS Reactor Physics Division Eugene Wigner Award, and the University of Chicago Distinguished Performance Award. In spite of his great stature and many accomplishments, Ely was a mild and modest gentleman who always gave full credit to others’ work, and was very approachable and an excellent listener. His technical questions at meetings were insightful, probing, and gentle. He also pursued the understanding of others’ points of view in personal and political matters with both intellect and sensitivity. His restaurant adventures at meetings and other venues have provided a rich array of gastronomic experiences and many fond memories to his many friends in our profession.

The Gelbard Review Lecture Series Conducted during the American Nuclear Society’s Conference Nuclear Mathematical and Computational Sciences: A Century in Review, A Century Anew Gatlinburg, Tennessee, April 6–10, 2003

Back row: Richard Hwang, Elmer Lewis, Kord Smith, Enrico Sartori Front row: Yousry Azmy, Elliott Whitesides, Jerry Spanier, Jack Dorning, Ed Larsen

Contents

Preface .. . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . vii Obituary Composed by Dr. Roger Blomquist for Dr. Ely Meyer Gelbard .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . xi 1

Advances in Discrete-Ordinates Methodology .. . . . . . . . . . . . . . . . . .. . . . . . . . . . . Edward W. Larsen and Jim E. Morel

1

2

Second-Order Neutron Transport Methods . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . 85 E.E. Lewis

3

Monte Carlo Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . .117 Jerome Spanier

4

Reactor Core Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . .167 Robert Roy

5

Resonance Theory in Reactor Applications . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . .217 R.N. Hwang

6

Sensitivity and Uncertainty Analysis of Models and Data.. . . . .. . . . . . . . . . .291 Dan Gabriel Cacuci

7

Criticality Safety Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . .355 G.E. Whitesides, R.M. Westfall, and C.M. Hopper

8

Nuclear Reactor Kinetics: 1934–1999 and Beyond . . . . . . . . . . . . . .. . . . . . . . . . .375 Jack Dorning

Index . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . .459

xv

Chapter 1

Advances in Discrete-Ordinates Methodology Edward W. Larsen and Jim E. Morel

1.1 Introduction In 1968, Bengt Carlson and Kaye Lathrop published a comprehensive review on the state of the art in discrete-ordinates .SN / calculations [10]. At that time, SN methodology existed primarily for reactor physics simulations. By today’s standards, those capabilities were limited, due to the less-developed theoretical state of SN methods and the slower and smaller computers that were then available. In this chapter, we review some of the major advances in SN methodology that have occurred since 1968. These advances, combined with the faster speeds and larger memories of today’s computers, enable today’s SN codes to simulate problems of much greater complexity, realism, and physical variety. Since 1968, several books and reviews on general numerical methods for SN simulations have been published [32, 46, 71], but none of these covers the advanced work done during the past 20 years. The specific purpose of this chapter is to describe how the field of SN calculations has matured through the lens of three important physical problems that can be simulated today but could not be realistically simulated in 1968. By discussing these problems and the methods developed to overcome their calculational difficulties, we hope to (i) show how dramatically the field of SN simulations of the transport equation has advanced and (ii) provide an introduction to the new algorithmic techniques that have enabled these advances. An outline of the remainder of this review follows. In Section 1.2, we briefly introduce the transport equation and discuss its basic temporal (implicit), energy (multigroup), directional .SN /, and spatial (finite-difference) discretizations, together with iterative solution procedures – as of 1968. The purpose of this section is to establish notation and set the stage for the later sections, which describe more recent developments. E.W. Larsen () Department of Nuclear Engineering and Radiological Sciences, University of Michigan, Ann Arbor, 48109-2104 Michigan, USA e-mail: [email protected] J.E. Morel Department of Nuclear Engineering, Texas A&M University, College Station, Texas, USA e-mail: [email protected] Y. Azmy and E. Sartori, Nuclear Computational Science: A Century in Review, c Springer Science+Business Media B.V. 2010 DOI 10.1007/978-90-481-3411-3 1,

1

2

E.W. Larsen and J.E. Morel

Section 1.3 discusses three important physical problems that could not be simulated in 1968 but can be realistically simulated today: thermal radiation transport, charged-particle transport, and oil-well logging tool design. In Section 1.4, we discuss advanced spatial discretizations (characteristic methods, discontinuous finite-element methods [DFEMs], and nodal methods) and the asymptotic thick diffusion limit (a technique to predict the validity of SN spatial discretizations for diffusive systems with optically thick spatial cells). Section 1.5 describes advances in discretizations of the angular derivatives associated with curvilinear geometries and treatments of anisotropic scattering. Section 1.6 covers advances in angular and energy discretizations for charged particles; Section 1.7 describes advances in time discretizations. In Section 1.8, we discuss major advances in iteration acceleration: diffusionsynthetic acceleration (DSA), linear multifrequency-grey acceleration for thermal radiation transport, fission source acceleration for time-dependent calculations, and upscatter acceleration. Section 1.9 outlines the recent application of preconditioned Krylov methods. Section 1.10 concludes with a brief discussion of challenges for the future: robust finite-element methods on nonorthogonal grids, positive and monotone methods, efficient parallel sweep algorithms for unstructured grids, further development of Krylov methods for solving the SN equations, methods for charged-particle calculations with pencil-beam sources, Galerkin quadrature with positive generalized weights, and ray-effect mitigation.

1.2 Basic Concepts The physical process discussed in this chapter is the interaction of radiation with matter (radiation transport, or particle transport). The archetypical equation that describes these interactions is the linear Boltzmann equation (LBE) [2, 3, 7, 13]: 1 @ .r; ; E; t / C r .r; ; E; t / C †t .r ; E/ .r ; ; E; t / v @t Z 1Z D †s r; 0 ; E 0 ! E .r; 0 ; E 0 ; t / d0 dE 0 0 4 Z Z .r; E/ 1 Q .r; E; t/ C †f r; E 0 : (1.1) r; 0 ; E 0 ; t d0 dE 0 C 4 4 0 4

In full generality, this equation has seven independent variables: three spatial variables .r/, two direction-of-flight (or angular) variables , energy .E/, and time .t/. Particle transport problems are difficult and costly to simulate because, in part, of the high dimensionality of phase space. In this section, we discuss the basic numerical methods used to solve Eq. (1.1) in the principal large computer codes of the 1960s [4, 5, 10, 14]. We assume that the reader understands the physical meaning and basic mathematical properties of each of the terms in Eq. (1.1), and we have used notation that is broadly standard. The discussion in this section is terse; we refer the reader to standard texts [13, 71] for details.

1

Advances in Discrete-Ordinates Methodology

3

The LBE given in Eq. (1.1) describes neutron transport with scattering and fission interactions. Variations of this equation primarily involve the types of interactions that are included. For instance, a gamma-ray transport equation would not have a fission term. Systems of coupled transport equations, each similar to Eq. (1.1), are required to describe the coupled transport of multiple types of particles, e.g., coupled neutron gamma-ray transport in which neutrons interact with nuclei to create gamma-rays and gamma-rays interact with nuclei to create neutrons. The principal computational difficulties associated with Eq. (1.1) are common to essentially all variations of this equation that are associated with different physical applications. In this chapter, we describe numerical methods in terms of Eq. (1.1), or simpler versions of that equation whenever possible. We consider variations of Eq. (1.1) that correspond to different physical applications only when necessary. To begin, we mention a few technical details. First, the differential scattering cross section is commonly written as a Legendre polynomial expansion: 1 X 2m C 1 Pm 0 †s;m r; E 0 ! E : †s r; ; E ! E D 4 mD0

0

0

(1.2)

The Legendre moments †s;m are typically calculated and stored for each material region. Also, initial and boundary conditions must be specified for Eq. (1.1). If V denotes the physical system and t D 0 is the initial time, then Eq. (1.1) holds for all r 2 V; 2 4; 0 < E < 1, and t > 0. At t D 0, must be fully specified in V : .r ; ; E; 0/ D

i

.r; ; E/; r 2 V; 2 4 ; 0 < E < 1:

(1.3)

Also, must be specified on the boundary @V for directions of flight pointing into V : .r ; ; E; t/ D

b

.r; ; E; t/; r 2 @V; n < 0; 0 < E < 1; 0 < t: (1.4)

Here, n is the unit outer normal vector at the boundary point r 2 @V . Many important algorithmic concepts can be explained most easily for problems with planar-geometry symmetry, in which the geometry and solution depend on only one spatial variable x and one angular variable D i . (The unit vector i points in the positive x-direction.) For a planar-geometry system 0 x X; Eq. (1.1) simplifies to 1 @ .x; ; E; t/ @ .x; ; E; t/ C C †t .x; E/ .x; ; E; t/ v @t @x Z 1Z 1 †s x; 0 ; ; E 0 ! E x; 0 ; E 0 ; t d0 dE 0 D 0

1

.x; E/ C 2

Z 0

1Z 1 1

†f x; E 0

Q.x; E; t/ x; 0 ; E 0 ; t d0 dE 0 C ; (1.5) 2

4


where 1 X 2m C 1 Pm ./Pm 0 †s;m x; E 0 ! E : (1.6) †s x; 0 ; ; E 0 ! E D 2 mD0

The initial condition for Eq. (1.5) is .x; ; E; 0/ D

i

.x; ; E/; 0 < x < X; 1 1; 0 < E < 1;

(1.7)

and the boundary conditions are .0; ; E; t/ D

l

.X; ; E; t/ D

.; E; t/; 0 < 1; 0 < E < 1; 0 < t; r

(1.8a)

.; E; t/; 1 < 0; 0 < E < 1; 0 t: (1.8b)

Because D vN , where v is the particle speed and N is the particle density, physically must be non-negative. If the cross sections, inhomogeneous source, initial conditions, and boundary conditions in Eqs. (1.5) through (1.8) are all non-negative (as they must be physically), then it can be shown that the solution of these equations is non-negative. However, the positivity of does not necessarily hold when approximations (discretizations) of the LBE are imposed. A desirable feature of a discretization for the LBE is that the resulting approximate solution should be positive – or nearly so. We now sketch the basic discretization and solution methods for Eqs. (1.5) through (1.8), which existed in computer codes in the late 1960s. We begin with the discretization of time. The most widely used time-discretization technique for transport problems, even today, is implicit time differencing. For a time interval tk1=2 0 (flow from left to right), Eq. (1.52) can be analytically solved for †t .xxj 1=2 /=n .x/ D e n

(1.52) n .x/:

† s †t .xxj 1=2 /=n C 1 e

j : n;j 1=2 2†t (1.53)

Using this expression, the flux exiting the j th cell is n;j C1=2

D e †t xj =n

n;j 1=2

† s C 1 e †t xj =n

j ; 2†t

(1.54)

and the cell-average flux is n;j

D

n 1 e †t xj =n n;j 1=2 †t xj † n s †t xj =n 1e C 1

j : †t xj 2†t

(1.55)

The exiting flux (Eq. (1.54)) is used as the incident flux for the next, i.e., .j C 1/, cell, and the cell-average flux (Eq. (1.55)) is folded into an array that, on completion of the transport sweep in all discrete directions, yields a new estimate for j .

1


25

The solution is considered to be converged when the new value of j differs from the old value by less than a prespecified convergence criterion for all j . This algorithm, with .x/ represented as a constant in each cell, is the StepCharacteristic (SC) method [12, 25]. For 1-D problems, the SC method can be formulated as a weighted-diamond scheme, i.e., of the form described by Eqs. (1.28) and (1.30). The SC method readily generalizes to multidimensional problems on irregular grids, but in these circumstances it cannot be formulated as a weighteddiamond scheme. The SC method is currently used in several multidimensional production neutron transport codes. In applications of these codes, the spatial grid is optically thin, and accurate solutions are obtained. However, the SC method is not accurate for problems demanding optically thick meshes, e.g., of the type described in Section 1.3. Thus, more complicated characteristic methods have been proposed and implemented, in which the scattering source is represented within each spatial cell as a linear, or even quadratic function of the spatial variables [25]. As can be expected, with each increase in the polynomial order of the representation of .x/: The accuracy of the resulting solution increases. The computational effort required to process the extra algebraic complexity

increases. The computer memory needed to store the extra problem unknowns increases.

1.4.2 Linear Discontinuous Method Among the most flexible and successful of the noncharacteristic methods are the discontinuous finite-element (DFE) methods [18, 70, 98, 105, 109, 116]. The linear discontinuous (LD) method is perhaps the archetypical method in this class. This method is based on representations of n .x/ and .x/ that are linear within each cell, but discontinuous at cell edges. In the LD method, Eq. (1.52) is approximated by n

d

n .x/

dx

C†t

n .x/ D

†s 2 .0/ .1/

j C .x xj / j ; xj 1=2 < x < xj C1=2 ; 2 xj (1.56)

where xj is the center of the j th cell; now the representation of .x/ requires two unknowns per cell, j.0/ and j.1/ . If the operator on the left side of Eq. (1.56) were inverted exactly, as described previously, we would obtain a (linear) characteristic method. However, discontinuous finite-element methods are based on an approximate, rather than an exact, inversion of this operator.

26


The LD method employs the following linear-discontinuous representation of n .x/:

n .x/

D

8 ˆ ˆ ˆ
0; xj 1=2 0 and from the right for n < 0, but is discontinuous otherwise. The unknowns nj and n;j ˙1=2 in Eq. (1.57) and .0/ .1/

j and j in Eq. (1.56) are related by

j.0/ D

N X

nj wn ;

(1.58a)

nD1 .1/

j

D

X

n;j C1=2

nj

wn C

n >0

To obtain equations for

X nj

n;j 1=2

wn : (1.58b)

n 0 obtain n xj

n;j C1=2

C

n;j 1=2

2

nj

C

†t 3

n;j C1=2

nj

D

†s .1/

: (1.61) 6 j

1


27

Equations ( 1.59) and ( 1.61) are two linear algebraic equations for nj and n;j ˙1=2 . [For n < 0, the term . n;j C1=2 nj / in Eq. (1.61) is replaced by . n;j n;j 1=2 /.] The resulting LD method can be generalized to multidimensional problems on structured (rectangular or orthogonal) and unstructured (nonorthogonal) spatial grids [97]. Variants of the LD method have also been developed, such as the lumped LD method, which is more robust for optically thick spatial cells but less accurate for optically thin cells, and various corner balance methods [76, 83, 90, 96]. In general, LD-like methods are much more accurate and robust for the difficult physical problems described in Section 1.3 than finite-difference methods. LD methods require greater computer arithmetic and storage than finite-difference methods, but their increased accuracy and robustness usually more than compensates for these disadvantages.

1.4.3 Nodal Methods Another class of spatial differencing techniques that have been developed and widely used in the nuclear reactor community are nodal methods [28–30, 43–45]. In essence, a nodal method approximates a multidimensional transport equation by a coupled system of 1-D transport equations. Discretization techniques that are highly accurate for 1-D can then be utilized. To illustrate, let us consider an .x; y/-geometry version of Eq. (1.51a) on a spatial cell (assuming for simplicity Q D 0): n

@

n .x; y/

C n

@

n .x; y/

C †t

n .x; y/

D

@x @y xi 1=2 < x < xi C1=2 ; yj 1=2 < y < yj C1=2 :

†s

.x; y/; 4 (1.62)

Transversely integrating this equation by the operator < >y;j D

1 yj

Z

yj C1=2

./ dy;

(1.63a)

>y;j ;

(1.63b)

yj 1=2

and defining n;y;j .x/

D
x;i D

1 xi

Z

xi C1=2

./dx;

(1.65a)

>x;i ;

(1.65b)

xi 1=2

and defining n;x;i .y/

D
0 and > 0 : n;y;j .x/

D

†t .x xi 1=2 / n †t .x xi 1=2 / 1 exp ; n

n;y;j .xi 1=2 / exp

C

qn;y;j;i †t

(1.70)

where qn;y;j;i D

†s 4 xi

Z

xi C1=2

y;j .x/dx xi 1=2

n yj

n;x;i .yj C1=2 /

n;x;i .yj 1=2 /

:

(1.71)

Following the same procedure for Eq. (1.69), we obtain †t .y yj 1=2 / n;x;i .y/ D n;x;i .yj 1=2 / exp n †t .y yj 1=2 / qn;x;i;j C 1 exp ; †t n

(1.72)

where qn;x;i;j D

†s 4 yj

Z

yj C1=2

x;i .y/ d y yj 1=2

n xi

n;y;j .xi C1=2 /

n;y;j .xi 1=2 /

:

(1.73)

Equation (1.70) is now evaluated at x D xi C1=2 , Eq. (1.72) is evaluated at y D yj C1=2 , and the resulting two equations [together with Eqs. (1.71) and (1.73)] are solved for the outgoing edge fluxes n;x;i .yj C1=2 / and n;y;j .xi C1=2 /. This procedure constitutes the 2-D Constant-Constant Nodal (CCN) method. This method is so-named because the transverse derivative term and the scattering source are approximated as constants in each transverse-integrated equation. The more accurate (but more expensive) 2-D Linear-Linear Nodal (LLN) method has four transverse-integrated equations with linear transverse derivatives and linear scattering sources in each equation. Extending these methods to 3-D is straightforward. Nodal transport (and diffusion) methods have played an extremely important role in the nuclear engineering community during the past 20 years. As discussed in Section 1.3.3 nodal methods (on rectangular cells) were applied to oil-well logging problems in the 1980s but were abandoned for that application because of their lack of suitability to nonorthogonal grids. The applicability of nodal methods to nonorthogonal grids remains a research topic.

30


1.4.4 Solution Accuracy in the Thick Diffusion Limit The last topic in this section is not a discretization scheme, but rather a theoretical technique, which, in the past 20 years, has become essential in predicting the accuracy of discretization schemes for diffusive problems with optically thick spatial cells .†t x 1/. Early (finite difference) spatial differencing schemes for the transport equation were experimentally and theoretically understood to be accurate only when spatial cells were optically thin .†t x 1/, and the accuracy of these schemes was generally measured by the order of their truncation error [25]. For example, an nth-order scheme would satisfy k

exact

x k

D O. n /; 1;

where jj jj is a suitable error norm and D †t x. In slab geometry, the DD and SC schemes are second-order, while the LD method is third-order. Analyses to mathematically prove the order of convergence were always carried out in slab geometry; in multidimensional geometries the SN solutions have singular characteristics, across which the solution is not smooth [20], so the truncation error analyses that can be carried out in slab geometry are not applicable. In fact, computer experiments have shown [33] that because of the singular characteristics, the order of convergence of the DD scheme in x,y-geometry depends on the definition of the error norm. Worse yet, the difficult thermal radiation and charged-particle transport problems described in Section 1.3 are so optically thick that it is impossible, because of limits in computer memory, to assign spatial grids for them that are optically thin. However, such calculations are associated with the diffusion and Fokker–Planck limits, respectively, and the spatial scale lengths for the solution associated with these limits are much larger than a mean free path. It is not unreasonable to expect a transport spatial discretization scheme to yield accurate results with optically thick cells if the scale length of the solution is well-resolved by the mesh. Indeed, one would intuitively expect to get accurate results with such mesh resolution. The difficulty is that a truncation error analysis does not provide useful information for these types of problems. Such an analysis tells us only that accurate results will be obtained by using optically thin cells. To determine if accurate results can be obtained with a mesh that is optically thick but resolves the spatial scale length of the solution, it is necessary to perform a discrete asymptotic analysis [49, 113]. Although the system is optically thick in both the diffusion and Fokker–Planck limits, the requirements associated with each limit for spatial differencing schemes are quite different. Accurate and robust spatial discretization schemes are generally required for charged-particle transport, but the highly anisotropic scattering treatment is of primary importance in the Fokker–Planck limit rather than the spatial differencing scheme [113]. In contrast, the spatial discretization is of primary importance in the diffusion limit. We focus on the diffusion limit here.

1


31

A diffusive problem is optically thick with weak absorption; it is a problem for which the transport solution is well-approximated by the diffusion solution. A spatial discretization of the SN equations is of practical use for diffusive problems if it possesses the optically thick diffusion limit [49, 52, 66, 86, 88, 105]. Such a discretization scheme will yield accurate results for diffusive problems if the spatial mesh cells are thin with respect to a diffusion length (the spatial scale length for the diffusion solution), even if these cells are thick with respect to a mean free path. To describe the diffusion limit, let us consider the monoenergetic planargeometry SN equations n

d

n .x/

dx

C †t

n .x/

D

N †s X 2 0

n0 .x/wn0 C

n D1

Q.x/ ; 1 n N; 2

(1.74)

and their diffusion approximation

d 1 d

.x/ C †a .x/ D Q.x/: dx 3†t dx

(1.75)

Here, we have used the standard notation †a D †t †s D absorption cross section; and

.x/ D

N X

n .x/wn

D scalar flux:

(1.76a)

(1.76b)

nD1

To motivate the subsequent analysis, we multiply the diffusion Eq. (1.75) by a positive constant ": d " d

.x/ C "†a .x/ D "Q.x/: (1.77) dx 3†t dx Clearly, the solution of the diffusion equation is unchanged. This shows that if we define the following scaled cross sections and source, †t ; " †a ! "†a ;

(1.78b)

Q.x/ ! "Q.x/;

(1.78c)

†t !

(1.78a)

which implies †s D †t †a !

†t "†a ; "

(1.78d)

32


then the diffusion equation is invariant under this scaling, for any choice of ". However, the SN equations are not invariant; they become n

d

n .x/

dx

C

†t "

1 2

n .x/ D

†t "†a "

Q.x/ ; 2

C"

X N

n0 .x/wn0

n0 D1

1 n N:

(1.79)

Now, one can show that for " 1, the solution of Eqs. (1.79) satisfies n .x/

.x/ C O."/; 2

D

(1.80)

where .x/ satisfies the diffusion Eq. (1.75). To derive this result, we solve Eqs. (1.79) by assuming a solution that, for " 1, depends on " by a simple asymptotic expansion: 1 X .x/ D "i n.i / .x/: (1.81) n i 0

Introducing Eq. (1.81) into Eq. (1.79) and equating the coefficients of different powers of ", we obtain for i 0 the following system of equations: †t

.i / n .x/

!

N 1 X 2 0

.i / 0 n0 .x/wn

n D1

D n .1/

d dx

.i 1/ .x/ n

N †a X 2 0

.i 2/ .x/ wn0 n0

n D1

C ıi;2

Q.x/ ; (1.82) 2

.2/

where n .x/ D n .x/ D 0. We solve this system recursively, by solving the first .i D 0/ equation, then the second .i D 1/ equation, etc. The first .i D 0/ equation is ! N 1 X .0/ .0/ 0 D 0: (1.83) †t n .x/ n0 .x/wn 2 0 n D1

Assuming that the quadrature set satisfies N X

Z wn D

nD1

1

d D 2;

(1.84)

1

Eq. (1.83) has the general isotropic solution .0/ n .x/

D

.0/ .x/ ; 2

where .0/ .x/ is – for now – undetermined.

(1.85)

1


33

The second .i D 1/ equation is, using Eq. (1.85), .1/ n .x/

†t

N 1 X 2 0

! .1/ 0 n0 .x/wn

D

n D1

n d .0/

.x/: 2 dx

(1.86)

Assuming that the quadrature set satisfies N X

Z n wn D

1

d D 0;

(1.87)

1

nD1

the general solution of Eq. (1.86) is .1/ n .x/

D

n d .0/

.1/ .x/ .x/; 2 2†t dx

(1.88)

where .1/ .x/ is undetermined. The third .i D 2/ equation, using Eqs. (1.85) and (1.86), is X t

.2/ n .x/

N 1 X 2 0

! .2/ 0 n0 .x/ wn

n D1

d D n dx

n d .0/ .x/

.1/ .x/ 2 2†t dx

!

Q.x/ †a .0/

.x/ C : 2 2

(1.89)

Unlike Eqs. (1.83) and (1.86), this third equation does not automatically possess a solution. To see this, we multiply Eq. (1.89) by wn and sum over 1 n N ; assuming that the quadrature set satisfies N X

Z 2n wn D

nD1

1

2 d D 1

2 ; 3

(1.90)

we obtain the solvability condition 0D

d 1 d .0/ .x/ †a .0/ .x/ C Q.x/: dx 3†t dx

(1.91)

If this equation is satisfied, then it can easily be shown that solutions to Eq. (1.89) exist. Equations (1.91), (1.85), and (1.81) confirm the result (1.80). The asymptotic analysis outlined above provides a direct mathematical link between the SN equations (1.74) and the diffusion equation (1.75). This analysis shows that if the cross sections and source in the SN equations are scaled by Eqs. (1.78) with " 1, the diffusion equation (1.75) is obtained. The condition " 1 is

34


consistent with the physical understanding of neutron diffusion: the mean free path [ D 1=†t D O."/] is small, the absorption rate [†a D O."/] is small, and the source [Q D O."/] is small – all of these “smallnesses” being balanced so that the resulting angular flux is O(1) and satisfies the diffusion equation. In the thick diffusion limit, two limits occur that have significance for spatial discretizations: 1 , 1. †t ! O "

.x/ . 2. n (x) ! 2 Thus, the total cross section becomes unbounded, yet the SN solution limits to an O.1/ diffusion solution. This result applies to the spatially continuous SN equations (no spatial discretization). We now ask: What happens if this same asymptotic analysis is applied to the spatially discretized SN equations? More precisely, let us consider a spatially discrete SN problem posed on a fixed spatial grid. We scale the cross sections and source in this problem exactly as in Eqs. (1.78). For " 1, we seek a solution of this discrete system in the form of Eqs. (1.81), i.e., we expand all unknowns (cell-average fluxes, cell-edge fluxes, etc.) as power series in ", and we solve the resulting hierarchy of equations as described above for the continuous SN equations. What happens to the spatially discrete SN solution in this limit? There are two possible answers to this question. First, because as " ! 0 the SN solution smoothly limits to the diffusion solution, it is plausible to hope that the spatially discrete SN solution will smoothly limit to the solution of a spatially discrete diffusion solution. (Then, if the chosen spatial grid is adequate to resolve the solution of this discrete diffusion problem, the resulting discrete solution will be accurate.) However, because †t D O."1 /, the optical thickness of spatial cells †t x ! 1 as " ! 0. This and the fact that SN solutions generally become inaccurate as †t x increases suggest that spatially discrete SN solutions may not limit to an accurate result as " ! 0. Which of these two possibilities is correct? The answer to this question depends on the chosen spatial discretization scheme. Some schemes are accurate in the thick diffusion limit; others are not. For example, the Step-Characteristic (SC) scheme fails as " ! 0 (the SC solution ! 0). The Diamond-Difference (DD) scheme fails unless all the diffusive regions of the problem have isotropic incident boundary fluxes (in the presence of nonisotropic boundary fluxes, DD solutions become corrupted by unphysical spatial oscillations). LD-like schemes perform successfully in the thick diffusion limit in 1-D geometries. LD methods also perform well in multi-D geometries with triangular (2-D) or tetrahedral (3-D) spatial grids, but they fail in quadrilateral (2-D) or hexahedral (3-D) grids. (However, bilinear-discontinuous methods work well for quadrilateral grids and trilinear-discontinuous methods work well for hexahedral grids.) The thick diffusion limit analysis, which has been applied to these discretization schemes and many others, accurately predicts the performance of approximation schemes in realistic calculations. This analysis has enabled the successful development of spatial discretization methods for problems with optically thick, diffusive

1


35

systems – in particular, for the thermal radiation transport and charged-particle transport problems discussed in Section 1.3. A reader may ask that if a transport problem is diffusive, then why not solve a simpler diffusion problem instead? The answer is that in many applications, only a part of the physical system is diffusive, and it may not be obvious where this diffusive part is. Also, some energy groups may be diffusive, while others are not. Finally, for time-dependent problems, some regions of space-energy phase space may be diffusive for certain times but not for others. For these reasons, it is generally infeasible to calculate accurate transport solutions by using the diffusion approximation in subregions of phase space where it is accurate. This leads to an important issue that can be discussed only briefly here: the behavior of SN spatial discretization schemes in the presence of unresolved boundary layers. (These are thin volumes, typically only a few mean free paths in width, containing the material boundaries that separate diffusive and nondiffusive subregions of a problem.) Across boundary layers, the flux usually has a rapid spatial variation; if the spatial grid is not sufficiently fine to resolve this fast variation, the boundary layer is said to be unresolved.) Many problems exist in which, due to computer memory limitations, it is not practical to prescribe a spatial grid that adequately resolves all boundary layers. Thus, one is led to the question of whether a given discretization scheme is accurate across an unresolved boundary layer. In particular, if an optically thick, diffusive region is adjacent to a nondiffusive region, can anything be said about the ability of a given discretization scheme to predict the changes in the flux across an unresolved boundary layer between two such regions? The asymptotic thick diffusion limit analysis does make it possible to study unresolved boundary layers; the conclusions so far are that no known differencing scheme is completely adequate to model unresolved boundary layers accurately. For example, LD methods are generally inaccurate in the first cell (containing the boundary layer) within the thick diffusive region, and they incorrectly predict that the flux exiting the diffusive region is isotropic. Generally, to be certain that a discrete solution is accurate, all spatial boundary layers must be adequately resolved by the spatial grid. For charged-particle transport problems, which are optically thick and have highly forward-peaked scattering, a more complicated asymptotic limit exists in which the total cross section †t D O."1 / and the mean scattering cosine 0 D 1 O."/. As " ! 0, the solution of the continuous transport equation limits to the solution of a Fokker–Planck equation (see the discussion in Section 1.3.2). Spaceangle discretization schemes have also been successfully analyzed in this asymptotic limit [113]. Ensuring that the discretized SN equations limit to a valid discretization of the Fokker–Planck equation is primarily related to the treatment of anisotropic scattering rather than the spatial differencing scheme. Nonetheless, the presence of very large and very small eigenvalues in the spectrum of the angular Fokker–Planck operator necessitates the use of accurate and robust spatial differencing schemes in Fokker–Planck calculations. In the thick diffusion and Fokker–Planck problems discussed above, it is generally impossible, given computer memory limitations, to use optically thin spatial

36


grids for the entire problem. To successfully simulate these problems, discretization schemes must produce accurate solutions for optically thick spatial grids away from boundary layers; and a theory is needed to justify the use of these schemes for these problems. The discontinuous finite-element schemes (such as LD and its variants) and the asymptotic (thick diffusion and Fokker–Planck) theories were developed to deal with just these practical difficulties.

1.5 Advances in Angular Discretizations Next we discuss (i) advances in SN discretizations for the angular derivative terms that appear in curvilinear coordinates systems and (ii) improvements to the standard SN treatment for highly anisotropic scattering. Perhaps surprisingly, very little has been accomplished during the past 40 years to successfully reduce the classic ray effects in SN simulations [21, 24].

1.5.1 Angular Derivatives The transport equation in curvilinear geometries contains one or more angular derivatives, in addition to spatial derivative terms. The traditional technique for treating the angular derivative term in the 1-D spherical geometry equation, which is representative of the traditional treatment used in essentially all curvilinear geometries, is described in Section 1.2. This technique is characterized by: The use of special ˇ-coefficients to represent the quantity .12 / at each angular

cell edge (see Eq. (1.23a)) The use of the diamond-in-angle relationship to express each cell-average angu-

lar flux in terms of the adjacent cell-edge angular fluxes (see Eq. (1.24)) The use of a starting-direction flux equation to obtain initial values for the angular

flux at D 1 (see Eq. (1.25))

This treatment has a deficiency, known as the discrete-ordinates flux dip, which consists of an erroneous suppression in the flux at the center of a sphere. Although the existence of the flux dip was recognized in the early 1960s, it was not eliminated until the early 1980s. Three features of the original method contributed to the existence of the flux dip: The starting-direction flux equation is a slab-geometry equation, but this was

originally put in the following curvilinear-like form before being spatially discretized [22]:

d 2 r dr D

1=2 N X n0 D1

C 2r

1=2

C †t .r/r 2

r 2 †s .r; 1; n0 /

1=2

n0 .r/wn0

C

r 2Q 2

(1.92)

1


37

This was thought to make the discretization for the starting-direction flux equation consistent with that of the other directions. However, it actually contributed to truncation errors that enhanced the flux dip. A boundary condition corresponding to specular reflection was used at the center

of the sphere, even though the angular flux at the center of a sphere is rigorously isotropic and equal to the starting-direction flux. This incorrectly allowed the angular flux at r D 0 to be anisotropic. The diamond-in-angle equation is inconsistent with the location of the quadrature cosines within each angular bin. As a result, the diamond-in-angle scheme does not preserve solutions that are linear in . In the late 1970s, it was proposed that the slab-geometry form of the startingdirection flux equation be discretized rather than the curvilinear-like form, and that all of the angular fluxes at the center of the sphere be set equal to the startingdirection flux value [26]. These two steps significantly reduced the severity of the flux dip. In the early 1980s, an angular weighted-diamond equation was proposed that related the angular edge and average fluxes in a manner consistent with the location of the cosine in each angular bin [38]: n .r/

D

n n1=2 wn

nC1=2 .r/

C

nC1=2 n wn

n1=2 .r/:

(1.93)

When all three of these measures were combined, the resulting angular discretization scheme eliminated the flux dip [38]. This scheme has been generalized to 2-D cylindrical geometry [38]. Very few practical improvements beyond the elimination of the flux dip have been made in SN angular derivative treatments. Discontinuous finite-element discretizations might have been expected to have had an impact, but this has not happened, partly because it is difficult to develop a discontinuous angular finiteelement method that is compatible with the standard SN method in multidimensional geometries. It is interesting that discontinuous angular derivative treatments do not require a starting-direction flux. This would appear to be an advantage, but one of the few linear-discontinuous SN angular derivative treatments ever developed for the 1-D spherical geometry equation was found to be less accurate than the weighteddiamond scheme (Eq. (1.93)) for a series of test problems [60]. A reason for this is that the starting-direction flux is computed (by Eq. (1.25)) with greater accuracy than the other directions; hence, significant accuracy is actually lost if the starting-direction flux plays no role in the angular derivative treatment. However, superior accuracy relative to the weighted-diamond scheme was obtained by using a quadratic-continuous approximation in the first angular cell and using a lineardiscontinuous approximation in the remaining angular cells [60]. All of these factors make it challenging to develop advanced SN angular derivative treatments [100].

38


1.5.2 Anisotropic Scattering The standard SN treatment for the scattering source, which is based on a Legendre polynomial expansion for the scattering cross section in conjunction with quadrature-generated spherical-harmonic moments of the angular flux, is still the workhorse for modern discrete-ordinates calculations, even though it is not always satisfactory. There are several reasons why it remains in widespread use: Fundamentally different approaches usually require significant processing of raw

cross-section data. Such techniques often have memory requirements that are significantly larger

than those of the standard treatment. The standard technique is often much more accurate than one would expect, even

when highly truncated cross-section expansions are used in a calculation. Next, we describe the standard method, together with an improvement that has had a notable impact on charged-particle calculations [54]. For simplicity, we consider the monoenergetic 1-D slab-geometry scattering source denoted by S : Z S .x; / D

C1 1

1 X 2m C 1 Pm ./Pm .0 /†s;m 2 mD0

!

x; 0 d0 :

(1.94)

We assume that the angular flux is a Legendre series of degree L, .x; / D where

L X 2m C 1 Pm ./ m .x/; 2 mD0

Z

m .x/ D

(1.95)

C1

1

Pm ./ .x; /d:

(1.96)

Substituting Eq. (1.95) into Eq. (1.94) and using the orthogonality of the Legendre polynomials, we find that the scattering source is S .x; / D

L X 2m C 1 Pm ./†s;m m .x/: 2 mD0

(1.97)

Thus, the scattering source generated by an angular flux that is a Legendre series of degree L is itself a Legendre series of degree no higher than L. Furthermore, the only cross-section information appearing in the scattering source is the first L C 1 moments of the scattering cross section. This same result is obtained if a cross-section expansion of degree L is used, rather than an exact expansion of infinite degree. In this case, the convergence of the cross-section expansion is irrelevant. This powerful result is not widely appreciated. An analogous result holds for multidimensional calculations when the angular flux takes the form of a finite

1


39

spherical-harmonic expansion. These results follow from the fact that the sphericalharmonic functions (which include the Legendre polynomials) are eigenfunctions of the Boltzmann scattering operator. We now discuss how this property impacts SN calculations. For simplicity, we consider the 1-D slab-geometry scattering source. Assuming that an N -point angular quadrature set is used in conjunction with a cross-section expansion of degree N 1, the SN scattering source takes the following form: S

n .x/ D

N 1 X

2m C 1 Pm .n /†s;m m .x/; 2 mD0

where

m .x/ D

N X

Pm .n /

n .x/wn :

(1.98)

(1.99)

nD1

We assume a Gauss–Legendre quadrature set. With N quadrature points, one can uniquely interpolate those points with a polynomial of degree N 1. Furthermore, since an N -point Gauss–Legendre set exactly integrates polynomials of degree 2N 1 [19], the Legendre moments in Eq. (1.99) are exactly the moments of the interpolatory polynomial. Considering our previous results regarding the scattering source for a polynomial angular flux representation, we see that the discrete scattering source values given in Eq. (1.98) are exactly those of the scattering source generated with the polynomial interpolation for the angular flux. Thus, if the true angular flux is well-represented by the polynomial interpolation of the discrete angular flux values, the true scattering source will similarly be well-represented by the polynomial interpolation of the discrete scattering source values. We again stress that this is true regardless of the convergence of the cross-section expansion. This property does not guarantee positive discrete scattering source values, given positive discrete angular flux values, because the polynomial interpolation of the discrete angular fluxes can be negative at some points, even though the discrete values themselves are positive. However, since polynomial interpolation at the Gauss points is known to be stable, any negativities in the angular flux interpolation will be small relative to the maximum discrete angular flux value. Therefore, any negativities in the discrete scattering source values will also be small relative to the maximum discrete scattering source value. Hence, accurate SN solutions for angle-integrated quantities can be obtained in a wide variety of problems in 1-D slab geometry with highly anisotropic scattering using Gauss–Legendre quadrature, even if the scattering cross-section expansion is highly truncated. If a Gauss–Legendre quadrature set is not used, some of the scattering source moments of the interpolatory polynomial will be properly computed, but others will not, depending on the accuracy of the quadrature set. It can be seen from Eq. (1.98) that the mth moment of the scattering source is just the product of the mth moment of the scattering cross section and the mth moment of the angular flux. Thus, any flux moment that is erroneous yields a corresponding scattering source moment that

40


is erroneous. This deficiency can be treated by generating a separate set of quadrature weights for each moment. In particular, for each 0 m N 1, one can generate N weights, fwm;n gN nD1 that are defined by the N linear equations N X

Pm .n /Pj .n /wm;n D

nD1

2 ımj ; 1 j N : 2m C 1

(1.100)

Pm .n /

(1.101)

Then Eq. (1.99) is replaced by

m .x/ D

N X

n .x/wn :

nD1

This method gives the desirable properties of Gauss quadrature to non-Gauss quadrature for the purpose of calculating the scattering source. (However, there is no guarantee that the weights generated in this way will be positive.) This is one variant of a more general technique known as Galerkin quadrature [54]. To present the more general method, we reexpress the standard SN technique for calculating the scattering source in terms of matrix algebra. In particular, we write Eqs. (1.98) and (1.99) as follows: SE D M†D E ;

(1.102)

where E is the vector of discrete angular flux values: E .

1;

2; : : : ;

N/

T

;

(1.103)

D is the N N matrix, Dm;n Pm .n / wn ;

(1.104)

† is the N N diagonal matrix: † diag.†0 ; †1 ; †2 ; : : :/;

(1.105)

and M is the N N matrix: Mn;m

2m C 1 Pm .n /: 2

(1.106)

The discrete-to-moment matrix D maps a vector of discrete angular flux values to a corresponding vector of Legendre flux moments. We note from Eq. (1.104) that the first row of this matrix consists of the standard quadrature weights, because P0 ./ D 1. The matrix † is the scattering matrix in the Legendre basis, or equivalently, the scattering matrix for the PN 1 approximation. It maps a vector of Legendre flux moments to a corresponding vector of Legendre scattering source

1


41

moments. The moment-to-discrete matrix M maps a vector of Legendre scattering source moments to a corresponding vector of discrete scattering source values. Using the orthogonal property of the Legendre polynomials, one can show that with Gauss quadrature, M D D1 , DM D

N X kD1

Di;k Mk;j D

N X kD1

Pi .k /

2j C 1 Pj .k /wk D ıi;j : 2

(1.107)

Thus, using Eq. (1.107), we can reexpress Eq. (1.102) as follows: SE D D1 †D E :

(1.108)

Equation (1.108) shows that the SN scattering matrix represents a similarity transformation of the Legendre scattering matrix, †. This means that the standard SN scattering source with Gauss quadrature (and a Legendre cross-section expansion of degree N 1) is equivalent to the scattering source of the PN 1 approximation. This is to be expected, considering the well-known equivalence between the SN and PN approximations in 1-D slab geometry [71]. If Gauss quadrature is not used, then M ¤ D1 , which is an undesirable result. The matrix D maps a vector of N discrete function values to N Legendre moments, and the matrix M maps a vector of N Legendre moments to N discrete function values. One can uniquely define a polynomial of degree N 1 either in terms of N Legendre moments or in terms of N discrete function values at N distinct points. Therefore, D and M should be inverses of one another. The moment-dependent weights defined in Eq. (1.100) ensure that this will be the case. We note that it is not necessary to actually generate the moment-dependent weights; one can directly obtain the correct matrix D simply by calculating the inverse of M. This Galerkin quadrature method is useful for 1-D calculations when quadratures with special directions are desired. For example, Lobatto and double Radau quadrature sets, which have quadrature points at D ˙1, are particularly useful for simulating a normally incident plane-wave of radiation [54]. In 2-D and 3-D, the Galerkin quadrature method is based on spherical-harmonic interpolation of the discrete angular fluxes rather than polynomial interpolation. Choosing the correct spherical harmonics for interpolation is more complicated in multidimensions because the number of spherical-harmonic functions of order N 1 does not equal the number of discrete directions in a multidimensional SN quadrature set. Nonetheless, suitable interpolation functions have been defined for triangular quadrature sets [54]. The Galerkin quadrature method in 2-D and 3-D can be much more accurate than the standard quadrature method with highly anisotropic scattering because there is no analog of Gauss quadrature in 2-D and 3-D, i.e., there is no 2-D or 3-D quadrature set that will exactly calculate all of the sphericalharmonic moments of the interpolated angular flux. In fact, fewer than half of the moments are exactly calculated with typical sets, e.g., even-moment symmetric sets [14], etc.

42


The Galerkin quadrature method can also accommodate nonpolynomial or nonspherical-harmonic interpolation functions [54]. To demonstrate this in 1-D, we consider a general interpolatory basis set for a given set of N discrete directions: N X

.x; / D

n .x/Bn ./;

(1.109)

nD1

where Bi .j / D ıij :

(1.110)

Multiplying Eq. (1.109) by Pm ./ and integrating over all directions, we obtain

m .x/ D

N X

Z n .x/

C1 1

nD1

Pm ./Bn ./d :

(1.111)

It follows from the definition of the discrete-to-moment matrix M and Eq. (1.110) that the components of D are Z Dm;n D

C1

1

Pm ./Bn ./d:

(1.112)

Equation (1.112) is valid for all types of interpolation functions, including polynomials. We note that the first row of the discrete-to-moment matrix consists of standard quadrature weights that are exact for integrating the interpolated angular flux: Z C1 N X .x; /d D (1.113)

.x/ D n wn ; 1

where

nD1

Z wn D

C1 1

Bn ./d:

(1.114)

These are called the companion quadrature weights. Nonpolynomial interpolation requires much more computational effort to generate the discrete-to-moment matrix, because the interpolatory basis functions must be explicitly formed and their products with the Legendre polynomials must be integrated. Also, one must invert the discrete-to-moment matrix to obtain the moment-to-discrete matrix, because the standard SN expression for the moment-to-discrete matrix, Eq. (1.106), is only correct for polynomial interpolation. We refer to the scattering source obtained by operating on the interpolated angular flux with the exact scattering kernel as the exact interpolation-generated scattering source. When the interpolation functions are nonpolynomial, the exact interpolation-generated scattering source is generally not expressible in terms of the interpolation functions. Thus, if the discrete scattering source values obtained from the Galerkin quadrature method are interpolated, one generally does not obtain

1


43

the exact interpolation-generated scattering source. Rather, one obtains a scattering source that has the same Legendre moments of degree 0 through N 1 as the exact interpolation-generated scattering source [54]. As an example of a useful nonpolynomial interpolation scheme, we consider a linear-discontinuous angular trial space in 1-D spherical geometry. Such a trial space is fully compatible with a linear-discontinuous treatment for the angular derivative term [60]. An “SN ” trial space of this type is defined to consist of N=2 equalwidth piecewise-linear segments in , where N is even and N > 2. There are two discrete angular flux unknowns per segment located p at the local Gauss S2 quadrature points, i.e., the points corresponding to ˙1= 3 obtained by linearly mapping Œ1; C1 onto each segment. A Galerkin quadrature set is generated for this trial space by exactly evaluating the Legendre angular flux moments of degree 0 through N 1 associated with the linear-discontinuous interpolation of the N discrete flux values [60]. The companion quadrature set corresponding to the Galerkin set, i.e., the standard quadrature set having the same quadrature points as the Galerkin set with quadrature weights that exactly integrate the interpolated angular flux representation, corresponds to a local Gauss S2 set on each linear segment. Since each local Gauss set exactly integrates cubic polynomials, it follows that the companion set will exactly evaluate the zeroth, first, and second Legendre moments of the interpolated angular flux. However, all higher flux moments will be inexactly evaluated, regardless of the quadrature order N . This is in contrast to the Galerkin quadrature set of order N , which always exactly evaluates the Legendre angular flux moments of degree 0 through N 1. Furthermore, because the companion quadrature set never exactly integrates polynomials of degree greater than 3, one cannot use a Legendre cross-section expansion of degree greater than 3 with the companion quadrature set (otherwise particle conservation will be lost). Thus, the accuracy of the scattering source with highly anisotropic scattering can be greatly improved for linear-discontinuous angular trial spaces in 1-D by using Galerkin quadrature. This enables one to use a linear-discontinuous approximation for the angular derivative term in 1-D spherical geometry in conjunction with an accurate treatment for highly anisotropic scattering. Perhaps the most important property of the Galerkin quadrature method, independent of the type of functions used to interpolate the discrete angular flux values, is that straight-ahead delta-function scattering is exactly treated. This has a very strong impact on charged-particle calculations, because it enables the total scattering cross section to be dramatically reduced (with an attendant decrease in the scattering ratio) while leaving the SN solution invariant. To demonstrate how straight-ahead scattering is exactly treated, let us consider the following differential scattering cross section (1.115) †s .0 / D ˛ı.0 1/; where ˛ is an arbitrary constant. The Boltzmann scattering operator associated with this cross section is ˛ times the identity operator: Z S

D

C1 1

˛ı.0 1/

0 0 d D ˛ ./:

(1.116)

44


Furthermore, the Legendre moments of this cross section are all equal to alpha: Z †m D ˛

C1 1

ı.0 1/Pm 0 d0 D ˛Pm .1/ D ˛; 0 m 1:

(1.117)

Thus, the diagonal matrix of cross-section moments used to construct the vector of discrete scattering source values is ˛ times the identity matrix: † D ˛I:

(1.118)

Substituting from Eq. (1.118) into Eq. (1.108), and recognizing that M D D1 , we obtain SE E D M˛ID E D ˛MD E D ˛ E ; (1.119) which agrees with Eq. (1.116). We have explicitly considered only the 1-D case, but this result also applies in multidimensions. For charged particles, the scattering ratio for each group is generally very close to unity, and the mean free path is very small; nonetheless, the transport process is not diffusive. This is because the “transport-corrected” scattering ratio, .†0 †1 /=†t , is not close to unity. Within-group straight-ahead scattering is equivalent to no scattering at all, since the particle scatters into the same group and direction it had before the scattering. Thus, one can add or subtract a straight-ahead differential scattering cross section from any physically correct within-group cross section without changing the analytic transport solution. Since all Galerkin quadratures treat straight-ahead scattering exactly, one can subtract the truncated expansion for a within-group straight-ahead scattering cross section from the physically correct cross-section expansion without changing the SN solution. For instance, let us consider the total Boltzmann scattering operator (outscatter minus inscatter) associated with a 1-D SN Galerkin quadrature: †0 E SE D ŒM†0 D M†D E D MŒ†s †D E D MŒdiag.†0 †0 ; †0 †1 ; : : : †0 †N 1 /D E : (1.120) Subtracting the delta-function cross section given in Eq. (1.115) from the physically correct cross section, we obtain the following modified outscatter matrix: †0 D .†0 ˛I/;

(1.121)

and the following modified inscatter matrix

where

.M†D/ D M† D;

(1.122)

† D diag.†0 ˛; †1 ˛; : : : ; †N 1 ˛/:

(1.123)

1


45

We note that while the outscatter and inscatter matrices are modified by subtraction of the straight-ahead scattering cross section, the total Boltzmann matrix, Eq. (1.120), does not change, i.e., †0 M† D D †0 M†D:

(1.124)

Thus, the SN solution does not change. However, the convergence properties of the source iteration process can be dramatically changed. A discussion of the optimal choice of ˛ is beyond the scope of this review, but the traditional choice (which is nearly optimal) is to set ˛ D †N . This extended transport correction can greatly reduce both the total scattering cross section and the scattering ratio in relativistic charged-particle transport calculations [9, 27]. We note that the significance of this cross-section modification depends on the convergence of the cross-section expansion. If the expansion is essentially converged, †N 1 will be very small relative to †0 , resulting in a negligible reduction in †0 ; and if the cross-section expansion is highly truncated, †N 1 will be comparable to †0 , resulting in a significant reduction in †0 . Acceptable computational efficiency often cannot be achieved without the use of the extended transport correction in charged-particle calculations. Thus, with Galerkin quadrature, the extended transport correction (correctly) leaves the SN solution invariant. This is a powerful motivation for using the Galerkin method in charged-particle calculations.

1.6 Advances in Fokker–Planck Discretizations Next, we discuss advances in Fokker–Planck angle and energy discretizations for charged-particle transport. We first consider the continuous-scattering operator, and then the continuous-slowing-down operator.

1.6.1 The Continuous-Scattering Operator The continuous-scattering operator in 1-D slab geometry is @ †r;t r @ 2 .x; ; E/: .x; ; E/ D 1 2 @ @

(1.125)

Taking the zeroth and first angular moments of Eq. (1.125), we obtain Z

C1 1

.x; / d D 0

(1.126)

46


and

Z

C1 1

.x; /d D †r;t r J.x/;

(1.127)

respectively, where J.x/ is the current Z J.x/ D

C1

.x; / d:

(1.128)

1

It is highly desirable for a numerical approximation to the continuous-scattering operator to preserve both the zeroth and first angular moments of that operator. Also, since the continuous-scattering operator is a diffusion operator on the unit sphere [see Section 1.3.2], it is highly desirable that the discretization for this operator yield a coefficient matrix that is symmetric and monotone. These two properties ensure that the matrix (like the analytic operator) will have only positive real eigenvalues and will yield positive solutions given positive sources. A straightforward discretization of Eq. (1.125) is †r;t r 1 nC1 n n n1 2 2 1 nC1=2 ; 1 n1=2 . /n D 2 wn nC1 n n n1 (1.129) where nC1=2 D n1=2 C wn ; 1 n N; 1=2 D 1:

(1.130)

This discretization results in a symmetric and monotone coefficient matrix and preserves Eq. (1.126) under numerical integration, but it preserves Eq. (1.127) only if each quadrature point lies at the center of its associated angular interval. (As previously noted, this never occurs with standard quadrature sets.) A discretization that does preserve Eq. (1.127) with standard quadrature sets is [40] †r;t r 1 . /n D 2 wn

ˇnC1=2

n n n1 ˇn1=2 nC1 n n n1 nC1

;

(1.131)

where the ˇ-coefficients are defined by Eq. (1.23b) and all else remains as previously defined. Thus the ˇ-coefficients used to admit the constant solution in the discretization of the angular derivative term in the spherical geometry transport equation are also used in the discretization of the continuous-scattering operator to preserve Eq. (1.127). This moment-preserving approach has been extended to multidimensions for product quadratures [120]. For instance, in three dimensions, the continuousscattering operator is

†r;t r D 2

@ 1 @2 2 @ .1 / C : @ @ 1 2 @! 2

(1.132)

1


47

Taking the zeroth and first angular moments of Eq. (1.132), we obtain Z

2

Z

C1

dd! D 0;

and

Z

2

Z

C1

dd! D †r;t r J ;

1

0

(1.133)

1

0

(1.134)

respectively, where J is the current Z

2

Z

C1

J D

0

dd!:

(1.135)

1

Standard triangular SN quadrature sets do not represent a rectangular angular mesh on the unit sphere, but product sets do. Hence, it is reasonably straightforward to derive a discretization for the continuous-scattering operator assuming a product quadrature set. An SN product quadrature set has 2N 2 directions and is formed by the tensor product of an N -point quadrature defined over the polar cosine and a 2N point quadrature defined over the azimuthal angle. Each direction can be uniquely referenced in terms of a polar index n and an azimuthal index j . A momentpreserving discretization for the 3-D continuous-scattering operator is 2 . /n;j D

†r;t r 6 4 2

1 p wn

ˇnC1=2

n 1 C 1 2 wa n j

nC1;j

ˇn1=2

n;j

nC1 n

n;j C1

n;j

!j C1 !j

n;j

n1;j

n n1

n;j

n;j 1

3 7 5:

!j !j 1

(1.136) p

Here, wn and waj are weights associated with the polar and azimuthal quadratures that sum to 2 and 2 , respectively, the ˇ-coefficients are identical to those defined for the 1-D case, and 2 K n ; 1 n 2N; 2N 1 cos N p Kn D 2.1 2n / C cn 1 2n; n D

ˇnC1=2 dnC1=2 ˇn1=2 dn1=2 ; p wn p p 1 2nC1 1 2n D : nC1 n

cn D dnC1=2

(1.137) (1.138) (1.139) (1.140)

48


The above discretization is defined only for product quadrature sets constructed with azimuthal quadrature sets of the Chebychev type: 2j 1 ; 1 j 2N; N ; 1 j 2N: waj D N !j D

(1.141) (1.142)

This discretization has a 5-point stencil, and is symmetric positive-definite (SPD) and monotone. It preserves Eqs. (1.133) and (1.134). The restriction to Chebychev azimuthal quadrature arises from the fact that three first-moment equations must be met, but the ˇ-coefficients and the coefficients can only be defined to meet two of them. The ˇ-coefficients are defined to preserve the moment equation associated with the polar cosine, i.e., cos , and for a general quadrature set, the -coefficients can preserve one of the two moment equations associated with the cosines that depend on the azimuthal angle, i.e., sin cos! or sin sin!. However, when a Chebychev azimuthal quadrature set is used, the -coefficients can be defined to preserve both azimuthal cosines. An alternative to a finite-difference representation for the continuous-scattering operator is a Legendre moment representation. Because the total Boltzmann scattering operator (outscatter minus inscatter) and the continuous-scattering operator have the same eigenfunctions, one can define effective cross-section moments to represent the continuous-scattering operator [31]. In particular, the mth eigenvalue of the continuous-scattering operator is c;m D †r;t r

m.m C 1/ ; 2

(1.143)

while the mth eigenvalue of the total Boltzmann scattering operator is b;m D †0 †m :

(1.144)

Without loss of generality, we assume that a 1-D SN calculation is performed with Galerkin quadrature. In this case, an effective cross-section expansion of degree N 1 is required. The first step in defining the effective cross-section moments is to equate the eigenvalues defined by Eqs. (1.143) and (1.144): †e0 †em D †r;t r

m.m C 1/ ; 0 m N 1: 2

(1.145)

This step does not uniquely define the effective cross-section moments, but rather leaves †0 a free parameter. The only consideration for choosing †0 is to minimize the effective scattering ratio. In analogy with the choice of ˛ in the extended transport correction,†0 is defined so that the last moment in the expansion is zero: †e0 D †r;t r

.N 1/N : 2

(1.146)

1


49

Substituting from Eq. (1.146) into Eq. (1.145), we obtain an expression for the remaining effective moments: †em D †r;t r

.N 1/N m.m C 1/ ; 1 m N 1: 2

(1.147)

When used in conjunction with Galerkin quadrature, the number of eigenvalues preserved is always equal to the number of discrete directions. This is to be contrasted with finite-difference approximations that only preserve the zeroth and first moments. Thus, the moment representation can be considered to be more accurate than the finite-difference representations, but the coefficient matrix for the moment representation is not monotone. Thus the moment representation is less robust than the finite-difference representation.

1.6.2 The Continuous-Slowing-Down Operator We next consider the continuous-slowing down operator: C .x; ; E/ D

@ˇ .x; ; E/ : @E

(1.148)

Discretizations of this operator have advanced from a multigroup or step-like treatment [31] through a diamond treatment [41] to a linear-discontinuous finite-element treatment [47]. Because standard SN codes were not originally intended to solve the charged-particle transport problems, early efforts in treating the continuousslowing-down operator were spent defining effective cross sections to implement various discretization schemes via the standard SN scattering source representation [31, 40, 47]. Modern codes that were designed to solve the charged-particle transport equation use a linear-discontinuous finite-element discretization in space, but treat the energy derivative similar to the spatial derivatives. Thus, this operator is inverted via a space-energy sweep [101]. Self-adjoint codes must still treat the term as a scattering source, and invert it via source iteration [112]. Because the spectral radius associated with iterations on the continuous-slowing down term can be very close to unity, these iterations must be accelerated. A synthetic acceleration technique, based on the diamond-difference approximation as the low-order operator, has been used for this purpose [47]. Although the pure Fokker–Planck equation can be solved, most charged-particle calculations are carried out with the Boltzmann–Fokker–Planck equation [37]. In this case, the Boltzmann scattering operator is treated with the standard multigroup approximation. The resulting hybrid multigroup/linear-discontinuous operator is formally treated as a linear-discontinuous operator. One simply takes the scattering kernel to be piecewise-constant in energy; equivalently, one takes all energy slopes associated with the scattering kernel to be zero.

50


For instance, let us assume that the angular flux within group g has the following linear-discontinuous dependence: .E/ D where

a;g

a;g

C

e;g

2 E Eg ; EgC1=2 E < Eg1=2 ; Eg

is the group average flux:

a;g

e;g

(1.149)

D

1 Eg

Z

Eg1=2

.E/ dE;

(1.150)

EgC1=2

is the group energy slope:

e;g

D

6 Eg2

Z

Eg1=2

E Eg .E/ dE;

(1.151)

EgC1=2

and Eg D Eg1=2 EgC1=2 is the group width for group g. We note from Eq. (1.149) that the angular flux at the interface energy between two groups is defined by the solution in the higher energy group. Since the continuous-slowing-down operator causes particles to lose energy, this choice is consistent with the direction of particle flow in energy. Under the assumptions of a Legendre expansion of degree L and 0 energy slopes for the multigroup scattering kernel, and a constant dependence of the restricted stopping power within each group, the discretized 1-D transport equation takes the following form for group g:

@

a;g

@x

C †t;g

a;g

D

G X L X 2m C 1 m Eg 0 †g 0 !g

a;m;g 0 Pm ./ C Qa;g 2 Eg 0 mD1

g D1

1 ˇr;g1 . a;g1 e;g1 / ˇr;g . a;g e;g / ; C Eg 3 @ e;g C †t;g e;g D Qe;g C ˇr;g1 a;g1 e;g1 @x Eg 2ˇr;g a;g C ˇr;g a;g

(1.152a)

e;g

:

(1.152b) Here, 'a;m;g denotes the group average Legendre flux moment of degree m for group g, Qa;g and Qe;g respectively denote the group average source and group source energy slope for group g, and †m g 0 !g is the standard multigroup-Legendre coefficient of degree m for a transfer from group g0 to group g: Z Eg1=2 Z Eg0 1=2 Z C1 1 D †s E 0 ! E; 0 Pm .0 /d0 dE 0 dE: †m g 0 !g Eg 0 EgC1=2 Eg0 C1=2 1 (1.153)

1


51

Equations (1.152a) and (1.152b) are solved via source iteration, but the continuousslowing-down operator is inverted during the sweep. In particular, the source iteration process takes the form

@

.`C1/ a;g

.`C1/ C †t;g a;g @x 1 h .`C1/ ˇr;g1 a;g1 Eg

D Qa;g C C @

.`C1/ e;g

@x

ˇr;g

.`C1/ a;g

.`C1/ e;g

i

L X 2m C 1 m .`/ †g!g a;m;g Pm ./ 2 mD1

G X g 0 DgC1

.`C1/ e;g1

C †t;g 2ˇr;g

L X 2m C 1 m Eg 0 .`C1/ †g 0 !g

0 Pm ./; 2 Eg a;m;g mD1

3 h ˇr;g1 Eg .`C1/ .`C1/ Cˇ r;g a;g a;g

.`C1/ e;g

.`C1/ a;g1

i

.`C1/ e;g

.`C1/ e;g1

(1.154a)

D Qe;g ; (1.154b)

where ` is the iteration index. If the full linear-discontinuous treatment for the scattering source were used, one would have nonzero scattering source energy slopes. In this case, one would iterate on both the scattering source averages and energy slopes. If the scattering ratio for a given group is sufficiently large to require convergence acceleration, the scattering source energy slopes could require acceleration in addition to the scattering source averages. Accelerating the source energy slopes significantly complicates the acceleration process. (We will address this point again, regarding the convergence acceleration of the temporal scattering source slopes associated with a linear-discontinuous discretization of the time derivative.) When combining the linear-discontinuous energy approximation with discontinuous finite-element approximations in space, one generally assumes a single energy slope per spatial cell rather than a separate energy slope for each spatial unknown within a cell. Although the latter assumption is more accurate, it can be excessively expensive. For example, if a trilinear-discontinuous spatial approximation is used for a 3-D rectangular cell, one gets eight spatial unknowns per cell per angle per group. If a separate energy slope is used for each spatial unknown, the number of unknowns per cell per angle per group increases to 16; but if only a single energy slope is used for the entire spatial cell, the number of unknowns only increases to nine. The linear-discontinuous discretization of the continuous-slowing-down operator represents a major improvement relative to step and diamond-difference discretizations [47]. In particular, the linear-discontinuous method is much less numerically diffusive than the step method and much less oscillatory than the diamond method. Nodal methods can also be applied to the continuous-slowing down operator. One would expect such methods to be comparable to discontinuous finite-element methods. However, because they have rarely been used in practice, we will not explicitly discuss nodal methods here.

52


1.7 Advances in Time Discretizations Next, we discuss advanced discretization techniques for the time derivative. As is the case for most of the derivatives terms in the Boltzmann equation, the time derivative has been treated with the discontinuous finite-element method [91,98] and the nodal method [59]. The linear-discontinuous method assumes an angular flux dependence of the following form over the kth time step: .t/ D where

k a

C

k a

k t

2 k t t ; t k1=2 < t t kC1=2 ; t k

is the average flux: k a

k t

(1.155)

D

1 t k

Z

t kC1=2

.t/dt ;

(1.156)

t tk .t/dt ;

(1.157)

t k1=2

is the temporal slope: k t

6 D .t k /2

Z

t kC1=2 t k1=2

t k D .t k1= 2 C t k1= 2 / = 2 is the midpoint of the time step, and t k D .t kC1= 2 t k1= 2 / is the width of the time step. We note from Eq. (1.155) that the angular flux at the interface between two time steps is defined by the solution from the previous time step. For simplicity, let us consider the 1-D time-dependent slab-geometry monoenergetic transport equation with isotropic scattering and an isotropic inhomogeneous source: @ Q †s 1@ C C †t D

C : (1.158) v @t @x 2 2 We obtain the following equations after applying the linear-discontinuous finiteelement approximation in time to Eq. (1.158): 1 h vt k

k a

C

k t

3 h vt k

k a

C

k t

D

†s k

C Qtk ; 2 t

2

k1 a

k a

C

C

k1 t

k1 a

i

C

C

k1 t

@ k C †t @x i

C

k

D

†s k

C Qak ; 2 (1.159a)

@ tk C †t;g @x

k t

(1.159b)

1


53

where Qak and Qtk respectively denote the source temporal average and source temporal slope for time step k. Equations (1.159a) and (1.159b) can be simultaneously solved via source iteration: i 1 h k;.`C1/ k;.`C1/ k1 k1 C C t a a t vt k k;.`C1/ @ a † s k;.`/ C †t ak;.`C1/ D

C C Qak ; (1.160a) @x 2 a i 3 h k;.`C1/ k;.`C1/ k;.`C1/ k1 k1 2 C C C t a a a t vt k k;.`C1/ @ †s k;.`/ C †t;g tk;.`C1/ D

C t C Qtk ; (1.160b) @x 2 t where ` is the iteration index. We note that one must iterate on both the temporal averages and the temporal slopes of the scattering source. If the scattering ratio is close to unity, the source iterations for both the averages and slopes must be accelerated. Deriving fully consistent diffusion acceleration equations from Eqs. (1.160a) and (1.160b) yields a complicated and difficult-to-solve system of coupled diffusion equations. If one uses step or diamond differencing in time, the diffusion-synthetic acceleration algorithm requires the solution of only one diffusion equation and is essentially identical to the algorithm for steady-state calculations. An approximate method has been developed in which the fully coupled system of diffusion acceleration equations associated with the linear-discontinuous temporal discretization scheme is replaced by two independent diffusion equations [98]. This approximate method appears to work quite well, resulting in a cost increase for performing DSA of about a factor of 2 relative to that associated with traditional temporal differencing schemes. While the linear-discontinuous finite-element approximation in time is more accurate than the step scheme and more robust than the diamond scheme, it is also more expensive. As with the continuous-slowing down operator, when one combines the linear-discontinuous temporal approximation with discontinuous finite-element approximations in space, one generally assumes a single temporal slope per spatial cell rather than a separate temporal slope for each spatial unknown within a cell. Although the latter assumption is more accurate, it can be excessively expensive. As we noted before, if a trilinear-discontinuous spatial approximation is used for a 3-D rectangular cell, one gets eight spatial unknowns per cell per angle per group. If a separate temporal slope is used for each spatial unknown, the number of unknowns per cell per angle per group increases to 16; but if only a single temporal slope for the entire spatial cell is used, the number of unknowns only increases to nine. In analogy with the derivation of Eqs. (1.70) and (1.72), we apply the constantconstant nodal method to Eq. (1.158) to obtain the following equations for > 0 (assuming for simplicity Q D 0): h i k1=2 k1=2 t exp † .t/ D v t t n;x;i n;x;i t h io qn;x;i;k n 1 exp †t v .t t k1=2 ; (1.161) C †t

54


where qn;x;i;k D

†s 4 t k

Z

t kC1=2 t k1=2

and n;t;k .x/ D

where †s qn;t;k;i D 4 xi

Z

n xi

n;t;k

xi C1=2

n;t;k .xi 1=2 /

t;k .x/ dx

;

(1.162)

†t .x xi 1=2 / .x / exp n;t;k x1=2 n †t .x xi 1=2 / qn;t;k;i C 1 exp ; †t n

xi C1=2 xi 1=2

x;i .t/ dt

1 h vt k

n;x;i

t kC1=2

n;x;i

(1.163)

i t k1=2 : (1.164)

The considerations for applying nodal methods in time are analogous to those for applying discontinuous finite-element methods. For instance, with a linear nodal method, one must be concerned with accelerating the temporal slopes of the scattering source. With a linear nodal method in both time and space [59], there would be only one temporal slope per space cell, but if one were to apply a nodal method in time in conjunction with another type of spatial discretization, multiple temporal slopes per space cell could arise. The practical need for advanced temporal discretization schemes relative to traditional discretization schemes is not as strong for the time variable as for the space and energy variables. This is due to the relative ease with which adaptive techniques can be applied to time integration, making it feasible to avoid the regimes in which simple discretization schemes perform poorly. In any event, the transport community has little experience with advanced temporal discretization schemes, and little research has been performed in this area.

1.8 Advances in Iteration Acceleration Next, we discuss major advances in iteration acceleration. A plethora of SN iterative acceleration techniques have been developed over the years (we refer the reader to the recent comprehensive review by Adams and Larsen [110]), but there is little doubt that the practical application of diffusion-synthetic acceleration (DSA) to source iteration has been the most significant advance in iteration acceleration techniques in the history of discrete-ordinates methods. Early attempts to utilize DSA were marred by successful performance only for problems with optically thin spatial grids [16]. Later, the subtleties concerning how these methods should be discretized became understood. To motivate the DSA and DSA-like methods, we first discuss the Fourier analysis technique, which has become an invaluable theoretical tool for predicting the convergence rate of iterative solutions of continuous and discrete problems.

1


55

1.8.1 Fourier Analysis The Source Iteration (SI) method is described in Section 1.2; see Eqs. (1.31) and (1.32). For a model infinite, homogeneous-medium transport problem with no discretization, 1 @ .x; / C †t .x; / D Œ†s .x/ C Q.x/; @x 2 Z 1

.x/ D x; 0 d0 ;

(1.165a) (1.165b)

1

the SI process begins with an initial guess .0/ .x/ of the scalar flux, and then for ` 1, the `th source iteration is defined by

@

.`/

.`1=2/

.x; /

C †t

@x .x/ D

.`1=2/

Z .`1=2/

.x; / D

1

.x/

.`1=2/

i 1h †s .`1/ .x/ C Q.x/ ; (1.166a) 2

x; 0 d0 :

(1.166b)

1

We now write an exact transport equation for the scalar and angular flux errors: ı .`1/ .x/ D .x/ .`1/ .x; /; ı

.`1=2/

.x; / D

.x; /

(1.167a)

.`1=2/

.x; /:

(1.167b)

To do this, we subtract Eqs. (1.166) from Eqs. (1.165), obtaining

@ı

ı

.`1=2/

.x; /

@x .`/

.x/ D ı

C †t ı

.`1=2/

Z .`1=2/

.x; / D

1

.x/

ı

.`1=2/

1 †s ı .`1/ .x/; (1.168a) 2

x; 0 d0 ;

(1.168b)

1

which define ı .`/ .x/ in terms of ı .`1/ .x/. Clearly, the rate of convergence of Eq. (1.166) is equal to the rate at which ı .`/ .x/ ! 0. To calculate this rate, we introduce the Fourier transforms Z 1 .`1/ ı

.x/ D a.`1/ ./e i †t x d; (1.169a) 1

ı

.`1=2/

.x; / D

Z

1 1

b .`1=2/ .; /e i †t x d;

(1.169b)

56


into Eqs. (1.168) to obtain c .i C 1/b .`1=2/ .; / D a.`1/ ./; 2 Z 1 a.`/ ./ D b .`1=2/ ; 0 d0 ;

(1.170a) (1.170b)

1

where c D †s =†t D scattering ratio. Equation (1.170a) gives b .`1=2/ .; / D

1 c a.`1/ ./ ; 2 1 C i

(1.171)

and then Eq. (1.170b) gives a.`/ ./ D

Z 1 d0 c a.`1/ ./ 2 1 1 C i 0

D ! ./ a.`1/ ./ D D Œ!./` a.0/ ./; where c !./ D 2

Z1 1

d0 c D tan1 ./ 1 C i 0

(1.172a)

(1.172b)

is the iteration eigenvalue. Equations (1.172) and (1.169a) yield Z1 ı

.`/

.x/ D

! ` ./a.0/ ./e i †t x d:

(1.173)

1

Thus, the rate at which the Fourier mode corresponding to wave number limits to zero is determined by !./. If j!./j 1, the corresponding mode converges rapidly. If j!./j < 1 and j!./j 1, the mode converges slowly. If j!./j 1, the mode does not converge. The overall rate of convergence is determined by the most slowly converging error mode, i.e., the largest value of j!./j over the Fourier variable . For large `, Eq. (1.173) implies

.`/

ı .x/ ` A;

(1.174)

where A is a constant, and

.`/

ı .x/

D sup j!./j D lim `!1 ı .`1/ .x/

1 0; n > 0; n > 0/, and assume that the cell-centered flux is constant inside the volume and that the surface (or cell-edged) fluxes are also constant over each of the cell’s six faces. The balance Eq. 4.21 can be approximated by n n n x ˆn C y ˆn C z ˆn D Qijk †ijk ˆn;ijk xi yj zk

(4.22)

where constant differences are taken for each coordinate: 8 C ˆ ˆ ˆ > ˆ > I < = X Q i .p/ K 1 ; P .u/ D hi .0/Gi .u/ Gi .u/ D L I ˆ > P ˆ > i D1 ˆ > Q hi .0/Ki .p/ ; :1

(5.73)

i D1

Hence, Eq. 5.72 provides the direct connection between the solution to the slowingdown equation without resonances to that with resonances. The Laplace transform of the scattering kernel can be determined analytically. The substitution of KO i .p/ into Eq. 5.71 yields, I P

LfP .u/g D

i D1

hi .0/ 1˛i

.p C 1/

I P i D1

1 e "i .1Cp/

hi .0/ 1˛i

1

e "i .1Cp/

(5.74)

By examining the above equation, the asymptotic properties of P .u/ and its shortrange nature of fluctuations can be readily established analytically. From the Laplace transform identity, it is quite obvious that the asymptotic limit of the Placzek function must be 1

lim P .u/ D lim Œp LfP .u/g D

u!0

p!0

1

I P i D1

D

˛i "i hi .0/ 1˛

1 N

(5.75)

i

where is the average lethargy increment per collision also defined in Eq. 5.67. Similarly, the short-range nature of fluctuations in P .u/ and/or Gi .u/ can also be determined via the use of the asymptotic expansion of Eq. 5.74 for large p. For illustration purposes, consider a single nuclide in the limit of large p, so that Eq. 5.74 can be represented by the following series, ˛ n n" p 1 e 1 1 ˛e "p X .1/n 1˛ lim LfP .u/g D (5.76) n ˛ ˛ p!1 1 ˛ p 1˛ nD0 p 1˛ where the subscript i is dropped for convenience. Note that the presence of exp.n"p/ indicates that its inverse must consist of a unit step function for each interval n". The inversion of Eq. 5.76 can be carried out term by term and, upon rearrangement, one obtains ˛ 1 e 1˛ u X .1/n ˛ n ˛ .un"/ C e 1˛ P .u/ D 1˛ nŠ 1 ˛ nD1

Œn C .u n"/ .u n"/n1 H.u n"/ where H.u n"/ is the Heaviside function (or unit step function).

(5.77)

246

R.N. Hwang

It is interesting to note that Eq. 5.77 provides the same results identifiable with those derived originally by Placzek if the sum in Eq. 5.77 is written out explicitly. A similar approach can be extended readily to P .u/ and/or Gi .u/ for a mixture with many nuclides. It can be shown that the results so obtained are identifiable with those given in [50]. Some results of Gi .u/ and P .u/ for a typical fast reactor composition are given both quantitatively as well as graphically in [50]. The fluctuations are noticeably less striking in the first few intervals when other scatterers with much smaller atomic weights are present. From Eq. 5.72, it is quite obvious that the fluctuations in the kernel Gi .uu0 / will impact the collision density as well when resonances are present. For the relatively high energy region where the extent of a resonance is usually less compared to the corresponding ", only the first scattering interval will be of practical interest from the perspective of resonance integral considerations. As a general rule, under the latter condition, F .u/ either exhibits a bump above its asymptotic value, or a drop below it, corresponding to whether gi .u/ assumes predominately positive or negative values. As discussed in [50], the former represents the scenario where the neutron width n is much greater than other partial widths whereas the latter represents the reverse scenario. The impact can be substantial when the neutron width is very large. For more details, the reader is referred to [50].

5.4.3 Resonance Integrals and Their Applications The resonance integral concept was introduced in the early stages of reactor physics development as an effective means to account for resonance absorption attributed to a handful of low-lying s-wave resonances of few actinides using the Breit–Wigner approximation. The general concept is readily extendable to other types of crosssection representations when cast into the form of the generalized pole expansion described in Section 5.2.

5.4.3.1 Traditional Resonance Integral Concept In addition to the use of the Breit–Wigner approximation, the three main assumptions used in earlier days that eventually led to four widely used approximations are: (1) flux recovery between resonances; (2) .0/ 1 above each resonance; (3) the collision density is taken to be constant for all nuclides in the mixture except for the resonant isotope in question. Hence, the slowing-down equation at energies far below the source energy becomes Zu †t .u/ .u/ D †m C u"

0

du0

e .uu / †sr .u0 / .u0 / 1˛

(5.78)

5

Resonance Theory in Reactor Applications

247

where †m is the energy-independent macroscopic scattering cross section of nonresonance nuclides and the subscripts for various parameters of the resonant isotope will, henceforth, be dropped for convenience. The simplified equation above also implies that the absorption of the medium is due entirely to the resonance under consideration. These assumptions lead immediately to various means to determine the flux and absorption rate for each Breit–Wigner resonance from which the slowing-down density can be determined in the context of Eq. 5.69. The resonance integral, defined as the absorption rate per atom, i.e., Z1 RI D

x .u/ .u/du

(5.79)

0

has been used as one of the major tools in studies of resonance absorption.

5.4.3.2 Various Resonance Integral Approximations During the period preceding the mid-1960s, there were four widely used approximations for treating the resonance integrals of the well-isolated Breit–Wigner resonances. Brief discussions of these methods are presented below.

Narrow Resonance Approximation (NR) If the extent of the resonance is small compared to the maximum energy loss per collision, the resonance contribution to the integral term in Eq. 5.78 becomes negligible so that Eq. 5.78 is reduced to

.u/ D

†p †p D †t .u/ †p C †Rt .u/

(5.80)

where †Rt .u/ is the total macroscopic resonance component and †p D †m C †pr is the total macroscopic potential scattering cross section. The corresponding resonance integral becomes p ar J.r ; ˇr ; ar / Er cos l Z1 .r ; x/dx p ar 1 D Er cos l 2 ˇr C .r ; x/ C ar .r ; x/

.RI/NR D

1

(5.81)

248

R.N. Hwang

where p D †p =Nr , ar D r C fr , ˇr D †p =.†0r cos l / and ar D tan 2 l . Here, †0r is the total macroscopic peak resonance cross section of level r. It is worth noting that the hard-sphere phase shift angle, l , was usually assumed to be small in the low energy range so that cos l 1 and sin 2 l 2 l . In fact, the asymmetric Doppler-broadened line-shape function was often ignored for the sake of expediency in many earlier works [6]. Since J.r ; ˇr /, in absence of ar , is readily amenable to the utilization of a precomputed table in a two-dimensional array, it was most widely used in earlier days. It will be shown that the integral of the form given in Eq. 5.81 can be computed efficiently without resorting to tabulation.

Infinite Mass Approximation (NRIM or WR) In the limit of infinite mass or the extent of the resonance much wider than the maximum lethargy increment per collision, i.e., " ! 0 in Eq. 5.78, the flux and its corresponding resonance integral are also reduced to the simple forms [6],

.u/ D

†m ; †m C †a .u/

.RI /NRIM D

m t J.r ; ˇr / Er

(5.82)

respectively, where ˇ r D †m =.†0r a =t / and m is the diluent cross section per absorber atom.

Intermediate Resonance Approximation (IR) It is apparent that some means to bridge the gap between the two approximations given above is needed. One widely used approximation that serves this purpose is the IR approximation originally proposed by Goldstein and Cohen [51]. From the foregoing discussions, it is reasonable to conjecture that the neutron flux across a resonance generally resembles the following approximate form,

.u/ D

†m C †pr †m C †a .u/ C †sr

(5.83)

where is a parameter characteristic of the resonance to be determined. This expression obviously leads to the NR and WR approximations as approaches 1 and 0 respectively. Furthermore, the corresponding resonance integral based on such a flux shape must retain the same general forms defined by Eqs. 5.81 and 5.82 provided that is insensitive to energy. By substituting the Doppler-broadened Breit–Wigner cross sections into Eq. 5.83, one obtains the same general form of the resonance integral as that for the NR-approximation with different arguments, i.e., .RI/IR D

./ p./ ar J r ; ˇr./ Er

(5.84)

5


249

./ where p./D.†m C †pr /=Nr , ar D ar t =.ar Cnr / and ˇr./ D ˇr t =.ar Cnr /, respectively. Thus, it amounts to redefining the parameters for the expression based on the NR-approximation defined by Eq. 5.81. Here, the parameter must reflect the higher-order effects of the slowing-down equation. Goldstein and Cohen [51] argued that Eq. 5.83 can be viewed as the first-order iterant of the integral equation of the Fredholm type defined by Eq. 5.78, i.e., .1/ D .u/. Substituting .1/ .u/ into Eq. 5.78, one obtains the second-order solution to be denoted by

.2/ .u/. If the iterative process converges rapidly, .1/ .u/ and .2/ .u/ must not be significantly different. Therefore, one plausible criterion for determining was by setting Z1 Z1 .1/ †ar .u/ .u/ du D †ar .u/ .2/ .u/ du (5.85) 0

0

From this transcendental equation, one may deduce the value of provided that the integration on the right-hand side can be carried out analytically into a manageable form. One obvious problem that hinders such a procedure is its complexity when the Doppler-broadened function along with its asymmetric component are considered. If one neglects the inherent temperature-dependence and uses the Lorentzian shape for the resonance, a procedure commonly used, can be determined explicitly.

Direct Numerical Approach (Nordheim’s Method) In the mid-1960s, the availability of modern computers made possible the alternative of a relatively rigorous treatment of the isolated resonance integral without resorting to significant approximations and complications illustrated by the IR approximation. One such method widely used in thermal reactor applications, especially in the United States, is that pioneered by Nordheim [52]. Nordheim’s method was intended for the treatment of the resonance integral in two-region repeated cells imbedded in an infinite reactor lattice via the collision probability method. For the purposes of this section, it suffices to focus only on its basic algorithm of treating the slowingdown equation and resonance integral for the case of the infinite homogeneous medium. The subject of collision probabilities will be addressed separately in the next section. This method amounts to a breakup of the resonance integral into the following form, Z1 Zu2 (5.86) .RI/Nord D x .u/ .u/du D x .u/ .u/du C I 0

u1

where the lethargies u1 and u2 correspond to the predetermined energy boundaries E1 and E2 , respectively, around the resonance peak with Ej D E0 ˙ m. p =2/; j 2 1; 2. Physically, the interval is taken to be an “m” multiple of the “practical” width p D t =2ˇ 1=2 approximately equivalent to the “half width” of the integrand based

250

R.N. Hwang

on the NR or WR approximations. If m p is taken to be much greater than the Doppler width, where the Doppler-broadened line-shape functions approach their Lorentzian limits, I , consisting of two integrals corresponding to those attributed to the tails on both sides of the resonance outside of u1 < u < u2 , becomes analytically integrable. To evaluate the integral defined in Eq. 5.86 requires the solution of the slowingdown equation given by Eq. 5.78 at each mesh point uj between u1 and u2 as well as the subsequent resonance integral itself. Thus, it amounts to the evaluation of a double integral which can be accomplished via the use of Simpson’s integration scheme with equally spaced mesh points with the mesh spacing taken to be much smaller than the maximum increment of lethargy per collision for the resonance absorber in question. Nordheim [52] showed that the discretized integrand for the integral defined in Eq. 5.86 can be readily obtained recursively. Unlike other approximations, this method preserves the rigor of the flux behavior within the crucial region around the peak of the Breit–Wigner resonance. It was eventually superseded by other more rigorous methods during the emergence of the fast reactor program when more advanced computational tools became available.

5.4.4 Various Developments Motivated by the Emergence of the Fast Reactor Program The emergence of fast reactor development and modern computing facilities had a significant impact on our philosophies for treating the resonance phenomena in practical applications. As discussed in Section 5.1.1, the former casts the role of resonance cross sections into a somewhat different light while the latter makes possible the development of many rigorous methods for routine applications unimaginable in earlier days. There were several challenges directly associated with fast reactor development. For resonance treatments in the homogeneous media, there are three major challenges. First, one must account for overlap effects resulting from neighboring resonances either from the same nuclide or from different nuclides even if the Breit– Wigner approximation is assumed. Second, one must also account for the spectral effects resulting from resonance scattering of intermediate weight nuclides. Third, one must deal with compatibility issues of the vastly improved nuclear data and their representations in conjunction with the existing reactor physics concepts and codes based on the Breit–Wigner approximation. As discussed in Section 5.2.3, the latter can be resolved via the use of the generalized pole representation. Hence, it suffices to focus on various developments in the first area for our purposes here. 5.4.4.1 Generalization and Computation of the J -Integral Our shift of interest to the higher energy region clearly justifies the extensive use of the NR-approximation for treating many sharp resonances of actinides. One

5


251

way to handle these issues is to generalize the integral within the context of the NR-approximation. If the total macroscopic resonance cross section defined in Eq. 5.80 is taken to be a linear combination of all Breit–Wigner resonances in the system, the NR-approximation at the k-th resonance can be cast into the following form via partial fractions [53], 0 B B 1

.u/ D †p B B † C † .u/ tk @ p P

1

C C !C C P A †p C †t k .u/ †p C †t k .u/ C †k 0 k .u/ k 0 ¤k

†k 0 k .u/

(5.87)

k 0 ¤k

where k 0 denotes the neighboring resonances. Thus, the corresponding generalized J -integral, referred to as the J -integral, must exhibit the following form, Jk D Jk .k ; ˇk ; ak ; bk /

X

Okk0

(5.88)

k 0 ¤k

where 1 Jk .k ; ˇk ; ak ; bk / D 2

Z1 1

.k ; xk / C bk .k ; xk / dxk ˇk C .k ; xk / C ak .k ; xk /

(5.89)

and Okk0 is referred to as the overlap term attributed to the k-th resonance with the integrand deducible from Eq. 5.87 in terms of the Doppler-broadened line-shape functions. The principal component J.k ; ˇk ; ak ; bk / represents the contribution from the k-th resonance alone. The parameter bk D 0 is taken when used in conjunction with the capture and fission resonance integrals so that it is consistent with Eq. 5.81. For total and/or scattering resonance integrals, one sets bk D ak . For details, the reader is referred to discussions given in [53]. It is important to point out that such an integral is composition-dependent and the required calculations must be carried out at run time. Hence, an efficient method for such a purpose is obviously needed. One proven method that has been extensively used in fast reactor applications is that based on the special form of the Gauss–Jacobi quadrature, sometimes also referred to as the Gauss–Chebyshev quadrature. While the detailed discussions are given in [53], some general features of this method will be discussed briefly for our purposes here. The rationale for this approach is based on the fact that the Doppler-broadened line-shape functions can become relatively well-behaved in a new domain upon

252

R.N. Hwang

changing the variable of integration. One simple way to accomplish this is via either one of the following two types of rational transformations [53], or u2 D C 2 xk2 = 1 C C 2 xk2 (5.90) u2 D 1= 1 C C 2 xk2 where C is a constant to be chosen. With no loss of generality, the integral of the following form can be cast into the form readily amenable to Gauss–Jacobi quadrature via the first type of transformation, i.e., Z1 f . .xk ; k /; .xk ; k //dxk 1

Z1 D 1

D

du p 1 u2

1 f . .xk ; k / ; .xk ; k / ; k / C u2

N 1 X f . .xk .un /; k /; .xk .un /; k /; k / C RN C N nD1 u2n

(5.91)

where un D cosŒ.2n 1/=2N and N is the total number of points considered. For such a quadrature to be effective, the part of the integrand inside the bracket in the equation above must be a smooth function of the new variable u so that it can be accurately approximated by a low-order Chebyshev polynomial expansion. For computation of the integral J.k ; ˇk ; ak ; bk /, it was shown in [53] that only a few mesh points are required to ensure accurate results if the constant C is carefully chosen in various ranges of k and ˇk . The same algorithm is equally applicable to the evaluation of the overlap correction term usually consisting of few contributing neighboring resonances if a larger number of mesh points is used. This algorithm has been incorporated into the MC2 -2 code [54] for routine applications. 5.4.4.2 Connections Between the Resonance Integral and Traditional Multigroup Cross Section Processing With no loss of generality, the multigroup cross section for a given reaction process x in a predetermined lethargy group ug D ug ug1 can be expressed in terms of resonance integrals within this group boundary as follows [53], P xk Jxk Fk =E0k .g/ p k2g .g/ Q x D (5.92) ug fg where fg , sometimes referred to as the flux correction factor, is defined as 1 X tk Jtk Fk fg D hFk i ug E0k

(5.93)

k2g

and Fk , the collision density at the k-th resonance, is generally taken to be constant.

5


253

There are two ways that the multigroup cross sections can be deployed in reactor applications. One simple approach known as the Bondarenko scheme where the selfshielding factors for a chosen group structure are precomputed by using Eq. 5.92, .g/ without the overlap term, as a function of p at various temperatures with results tabulated and stored in one-dimensional arrays. In lieu of the overlap term, .g/ the quantity p infers the approximate composition dependence of the multigroup cross sections so that this scheme provides an efficient means for survey-type calculations. Another approach is to process these multigroup cross sections taking into account the actual composition dependence before the results are passed on to neutronic calculations [54]. As mentioned earlier, accurate accounts of the composition involving many nuclides with resonance structures are extremely important in fast reactor applications and the multigroup cross sections must undergo further refinements prior to their deployment. For our purposes here, it suffices to present the rationale for two widely used methods that provide a better account of the global spectral effects.

Ultra-Fine Group/Fundamental Mode Approach The ultra-fine group (ufg) approach was pioneered by Hummel [12] in order to account for the fine spectrum effects attributed to resonances of nuclides with intermediate atomic weights not accounted for in the earlier resonance theory. It was intended specifically for the treatment of resonance structures of the structural and metal coolant materials present in a fast reactor composition. The dominant resonances of these materials are characterized by relatively wide neutron widths not sensitive to Doppler-broadening. One conceptually simple means to deal with this issue is to divide the entire energy span of interest into 2,000 ufgs with equally spaced lethargy widths u D 1=120, approximately equivalent to half of the maximum lethargy increment per collision of an actinide. One way to account for fine structure effects on the broad group cross sections is to use the ufg flux computed based on the predetermined ufg cross sections as the tool to collapse the latter into the desired group structure. This can be accomplished in two steps. First, compute the ufg cross sections according to Eq. 5.92. Second, once the ufg cross section set for a given composition is generated, the ufg spectrum can be computed via the usual fundamental mode spectrum calculations utilizing the consistent or inconsistent PN or BN approximation for diffusion equations widely used in reactor applications. This algorithm was coded in the MC2 -2 code [54].

Continuous Slowing-Down Approach The continuous slowing-down approach is an alternative in which the ufg weighting spectrum is determined by solving the slowing-down density equation given by Eq. 5.66. It is based on the rationale that the effects attributed to the relatively

254

R.N. Hwang

smooth-varying cross sections and those attributed to the sharp resonances can be treated separately, a method particularly amenable in conjunction with the resonance integral concept. If the slowing-down density defined by Eq. 5.66 can be determined in the absence of sharp resonances, the corresponding local slowing-down density with sharp resonances and thus the local flux can also be specified via the attenuation of the former using the resonance escape probability in the same context defined by Eq. 5.69. A comprehensive review of this subject was presented by Stacey [55]. For our purposes here, it suffices to focus only on the conceptual basis of this approach. One best known approximation for solving Eq. 5.66 for the case of relatively slow-varying cross sections was pioneered by Goertzel and Greuling [56] using the synthetic kernel approach. Their rationale can also be viewed as a natural consequence of applying a low-order Taylor’s expansion to the quantity †s .u/ .u/, the scattering component of the collision density [55]. The substitution of the expansion into Eq. 5.66 yields ! 1 X X .1/n .i / d n mn Œ†si .u/ .u/ q.u/ D nŠ dun nD0 i

Zu / m.i n D

.u u0 /n ki .u u0 /du0

(5.94)

u"i / where ki .u u0 / is defined in Eq. 5.66 and m.i n is the n-th order moment for the i-th nuclide. By retaining the first two terms in the expansion and the relation given by Eq. 5.68, the resulting first-order differential equation of q.u/ can be represented in .i / terms of moderating parameters that depend on the low-order mn and local cross sections. In absence of sharp resonances, such an equation can be solved readily. Because of its importance in fast reactor applications, a great deal of improvement of the original version by Goertzel and Greuling [56] has since been added, most notably by Stacey [55]. The improved version of Stacey [55] was based on a somewhat different low-order Taylor’s expansion from that given by Eq. 5.94. Instead of †s .u/ .u/, a similar expansion was made on †t .u/ .u/, the total collision density, which conceptually exhibits smoother behavior than the former. By so doing, the same type of first-order differential equation as the former was derived except that all moderating parameters must be redefined. For fast reactor applications, higher-order Legendre moments were also included in computing the moderating parameters. In particular, this improved method has been incorporated into the MC2 -2 code [54] as an option for computing the ufg fundamental mode spectrum in the resolved energy range. For details, the reader is referred to [54] and [55].

5.4.4.3 Rigorous Treatment of Resonance Absorption via Numerical Means There exist various degrees of inherent limitations in all resonance integral methods described so far. In general, there are two limitations in common. First, the rigor in the treatment of the slowing-down equation is lacking especially when resonances

5


255

of many nuclides are present. Second, the range of integration is generally taken from minus infinity to infinity without due consideration of the finite group structures in conjunction with the multigroup applications downstream. This gives rise to the so-called boundary effects. These limitations along with the need for accurate treatment of heterogeneous effects of reactor lattices to be addressed in the next section provided a strong motivation for the development of more rigorous methods particularly useful for benchmarking purposes. For our purpose here, one rigorous numerical treatment of the resonance absorption in homogeneous media will be presented. The method pioneered by Kier [57] and later improved by Olson [58], sometimes referred to as the “hyper-fine group” method, is believed to be tailor-made for such a purpose. For illustration purposes, first let us consider the simplest case involving only one nuclide. The rationale is to divide a given scattering interval into many equally spaced hyper-fine groups (hfg’s) with spacing u much smaller than the extent of the resonances under consideration. This allows us to describe the flux defined by Eq. 5.63 for any given group, say k , discretely if the elastic scattering process for the hfgs within a span of "i is known. Thus, it suffices to isolate the neutronic balance for one scattering interval in the presence of a single nuclide for illustration purposes. Let the total number of hfgs within the span of " below the hfg in question be L D "=u;

˛ D exp.Lu/

(5.95)

The neutronic balance in the hfg sense can be achieved via the use of the “effective” scattering kernel according to Kier [57]. Physically, if the scattering kernel K.u u0 / is viewed as the probability density function (p.d.f.), Kier’s “effective” scattering kernel [57] can be construed as the corresponding cumulative density function (c.d.f.) as a function of the successive hfgs within " above the hfg in question. If Pl denotes the c.d.f. for the l-th hfg above the initial group with lower boundary u0 , one can express it as Pl u D KQ l u; D KQ LST u C KQ s u;

1 l L1 lDL

(5.96)

where the c.d.f. for l D L represents the probability of neutrons initiated at u0 reaching the hfg in question and must take into account in-group scattering of that hfg in order to preserve neutron balance. These effective scattering kernels can be specified in the form of double integrals,

1 KQ l u D 1˛

u0ZCu

.l1/

Z

u0

u

du u0

ul0 u

0 e .uu / du0 D KQ 1 u e .li/u

(5.97)

256

R.N. Hwang

KQ LST u D

1 1˛

u0ZCu

u0 .L1/u Z

u0

KQ s u D

1 1˛

0 e .uu / du0 D KQ L u ˛ KQ s u (5.98)

du uLu

uLZCu

Zu

du uL

1 .u 1 C e u / 1˛

0

e .uu / du0 D

(5.99)

u0

Here, KQ s u is the probability the scattering is in-group. It can be readily shown that the sum of Pl over all l is equal to unity as expected. If one denotes the flux per unit lethargy for a given group k by k and the corresponding collision density per unit lethargy by Sk , the corresponding slowing-down equation can be discretized as follows. Sk D

L X

! Q Kl .†s /kl ˛Ps .†s /kL C .†s /k KQ s D Sk.in/ C Sk.s/ (5.100)

lD1

where Sk.in/ and Sk.s/ physically signify the in-scattering source from other hfgs and the self-scattering source of group k, respectively. One potential issue is the evaluation of the in-scattering source when an exceedingly large number of hfgs are considered. Kier [57] showed that the computational efficiency could be significantly enhanced via the use of the following recurrence relation derivable from the properties of KQ l given by Eq. 5.97, i.e., .in/ Sk.in/ D Sk1 .KQ 1 KQ s e u /˛.†s /k1L CKQ 1 .†s /k1 ˛ KQ s .†s /kL

(5.101)

Hence, for a scattering interval, one only needs to deal with four terms coming from the lethargy groups below k (or energy groups above k). For practical applications in conjunction with the multigroup approach, the calculation usually begins in the energy region slightly above the resolved energy region where the cross sections are taken to be constant. .CR/k D

Sk.in/

1 .†sk =†tk /KQ s

;

k D

.CR/k †tk

(5.102)

Once Sk.in/ is known, the corresponding collision rate, .CR/k , and flux, k , for the given k-th hfg can also be defined. The procedure described above can be repeated until the hfg fluxes and the corresponding reaction rates are determined. The effective group cross section for a lethargy group, u1 u2 , consisting of N hyper-fine groups is simply, !, N ! N X X (5.103) Q x D xk k

k nD1

nD1

5


257

The same procedure can be readily generalized for the practical case involving many nuclides. This approach is, in principle, rigorous so long as the slowing-down equation is valid. These basic algorithms have been incorporated into the RABBLEcode [57] and MC2 -2 code [54].

5.5 Resonance Absorption in Heterogeneous Media The treatment of resonance absorption in heterogeneous media obviously presents a greater challenge than that in homogeneous media. Unlike the latter, the energy and space dependence become intertwined. It is clearly unrealistic if the detailed distribution of neutron flux throughout the entire reactor must be specified simultaneously at any arbitrary energy and temperature. Hence, the general starting point for all existing deterministic approaches is to focus on the examination of the resonance absorption within a localized unit cell. A realistic reactor can be generally viewed as an ensemble of lattices consisting of unit cells. These unit cells, in principle, can have either identical or different composition depending on the design under consideration.

5.5.1 Traditional Collision Probability Methods for a Two-Region Cell One convenient way of treating resonance absorption in a reactor lattice is via the use of the collision probability method. The rationale of using such a method can be best illustrated via the two-region cell embedded in an infinite reactor lattice considered by Chernick [4, 5]. The cell consists of a fuel region and a large moderator region designated by 0 and 1, respectively and no inter-cell interactions are assumed. If F0 .u/ and F1 .u/ denote the collision density in these regions, respectively, the neutron conservation of this cell can be specified by the coupled integral equations given below, V0 F0 .u/ D .1 P0 /S0 C P1 S1

(5.104)

V1 F1 .u/ D .1 P1 /S1 C P0 S0

(5.105)

where the respective volumes and escape probabilities are denoted by V0 , V1 , P0 , and P1 , respectively. Alternatively, .1 Pn / is generally referred to as the collision probability for the n-th region. The quantity Sn represents the scattering source in the specific region and is given by Sn D

I.n/ X

Zu

i D1u"

i

.n/

Ki .u u0 /

†si

†.n/ t .u/

Vn Fn .u/d u0 ;

n 2 0 or 1

(5.106)

where I.n/ is the total number of nuclides in the nth region under consideration.

258

R.N. Hwang

Thus, the coupled equations amount to the combination of two slowing-down equations previously defined for homogeneous media connected together via the respective escape (or collision) probabilities. To specify the collision density as a function of u requires the explicit description of Pn as a function of lethargy (or energy) and the means to solve the coupled integral equations. In earlier days, the direct numerical solution to these equations, in conjunction with resonance integral calculations, was obviously prohibitive and the task required a great deal of simplifications. One general identity which is crucial to lattice physics studies is the widely used reciprocity relation. For a two-region cell, it is simply, .0/

.1/

P0 †t .u/ V0 D P1 †t .u/ V1

(5.107)

so that P1 can be computed once P0 is known. Another simplification is made possible by the fact that region 1 is usually assigned to the moderator with constant cross sections. Hence, it is reasonable to assume the NR-approximation in that region. By taking F1 .u/ D †1 .u/ and utilizing the reciprocity relation, the coupled integral equations can be reduced to a familiar form similar to that for homogeneous media described in the previous section, F0 .u/ D .1 P0 /

fuel X i

1 1 ˛i

Zu u"

0

e .uu /

0 †.0/ si .u / 0 †.0/ t .u /

F0 .u0 /d u0 CP0 †t .u/ (5.108) .0/

The simplified equation above is generally considered as the starting point for various approximations of the two-region cell configurations to follow. The key issue here is how to define the escape probability P0 as a function of †.0/ t .u/ so that the traditional approximations for homogeneous media can be applied consistently.

5.5.1.1 General Features of Collision Probability of Practical Interest For the sake of completeness, the transport theory origin of the escape probability method [59] will be briefly addressed. With no loss of generality, let us start with the generic representation of the escape probability for an arbitrary region with volume V and surface area S from which all subsequent simplifications are based, 0 1 Z Z h i 1 @ 1 E nE 0 / 1 exp.†t jErs rEs 0 j/d E AD 1T dErs 0 . Pesc D N S †t l †t lN S

(5.109) E where lN D 4V =S is commonly referred to as the average of the chord length, 0 is the unit vector along the direction of motion of the neutron, and nE is the unit vector normal to the surface S . The quantity inside the parentheses is denoted by

5


259

1 T , where T is known as the transmission probability signifying the fraction of neutrons transmitted through the region in question without suffering any collision. Alternatively, T also represents the neutron current of unit strength at the surface that passes through the region without collision. The corresponding collision probability is simply, (5.110) Pc D 1 Pesc Let l D jErs rEs 0 j be the chord length along the neutron path, a straight line connecting the intersections of the two vectors rEs and rEs 0 with the surface. For practical configurations of convex geometries, the outer integral over the surface in Eq. 5.109 becomes independent of the inner integral so that, Pesc D

S 4V †t .u/

Z

E nE / Œ1 exp.†t .u/l.d; // d .

(5.111)

where l.d; /, is a function of “diameter” (or thickness where appropriate) and the direction of flight . For the sake of clarity, the general attributes of geometry-related quantities in the integrand of Eq. 5.111 are given in Table 5.2 for the three most commonly used configurations in reactor physics applications. These are the fundamental geometric relationships in terms of ', the angle between the projection of l with respect to the x-axis, and , the azimuthal angle, used in defining their respective escape and/or collision probabilities. Let us begin with a simple scenario of an isolated unit cell in the sense that all neutrons escaping from the fuel lump will suffer their next collision in the surrounding moderator. For the isolated fuel configurations of practical interest, the escape probabilities are analytically derivable as illustrated in [59]. Figure 5.1 ilN a quantity commonly referred to lustrates their behavior as a function of †t .u/l, as the “optical” thickness, for slab, cylinder, and sphere unit cells. In spite of their differences in shape, the escape probabilities very much follow the same pattern. N In particular, they share ı the same limiting values, i.e., for †t l ! 0, Pesc D 1 and †t lN ! 1, Pesc D 1 †t lN known as the “white” and “black” limits, respectively. These limiting properties lead immediately to two well-known approximations, referred to as the mean-chord approximation and Wigner’s rational approximation [2, 3]. The former is defined as h i .†t .u/l/ (5.112) Pesc D 1 exp.†t .u/l/ Table 5.2 Geometry-dependent quantities required for computing collision probabilities E n E E Configurations .; E/ d l.d; / Infinite slab, d D t cos ' sin ' d' d t =sin ' Infinite rod, diam. d sin ' cos ' sin ' d' d d cos =sin ' Sphere, diameter ds cos ' sin ' d' d ds cos '

260

R.N. Hwang THREE CONFIGURATIONS & WIGNER APPROX.

ESCAPE PROBABILTY

1 cylinder slab sphere rational

0.8

0.6

0.4

0.2

0 2

0

4 6 8 OPTICAL THICKNESS

10

12

Fig. 5.1 Analytically based results and Wigner’s rational approximation

while the latter is defined as Pesc D 1

1 C †t .u/ l

(5.113)

It is quite apparent that these functions exhibit the same limiting properties as those of the rigorous representations. The discrepancies of the rational approximation, for instance, are usually no more than a few percent for the isolated lumps considered as illustrated also in Fig. 5.1. These approximations made possible much of the earlier work on lattice physics. 5.5.1.2 Various Earlier Methods Based on Approximate Escape Probabilities for Isolated Fuel Lumps The similarity between the simplified version of the slowing-down equation defined by Eq. 5.108 and that of the homogeneous media previously given makes possible the application of the comparable approximations for reactor physics calculations. For our purpose here, the most widely used NR and WR approximations will be presented for illustration purposes.

NR-Approximation .0/

Let F0 .u/ D †t .u/ 0 .u/ be the collision density for the fuel region. Under the NR-approximation as before, Eq. 5.108 for a single resonance absorber mixed with diluent becomes .0/ .0/ †.0/ t .u/ 0 .u/ D .1 P0 /†p C P0 †t .u/

(5.114)

5


261

Thus, the corresponding resonance integral can be expressed in terms of two integrals Z1 RI D

†.0/ p 0

a .u/ †.0/ t .u/

Z1 du C 0

P0 a.0/ .u/†.0/ Rt .u/ †.0/ t .u/

du D Iv C Is

(5.115)

.0/ where †.0/ p and †Rt .u/ are the total macroscopic potential scattering cross section and total macroscopic resonance cross section, respectively. The quantities Iv and Is are referred to as the “volume” integral and “surface” integral respectively according to Wigner et al. [3]. The former is identical to that for homogeneous media described in the previous section while the latter involves the escape probability characteristic of the heterogeneous nature of the cell. The general form of the surface integral becomes + *Z1 h i a.0/ .u/†.0/ .u/ S0 .0/ Rt 1 exp.†t .u/l/ du Is D (5.116) 4V0 †.0/ .u/ t 0

where the angular bracket denotes the integration over the geometric configuration in the context of Eq. 5.111. From the perspective of reactor applications, this approach is not widely used. In contrast, the simplest possible means of computing the resonance integral for the fuel lump is via the direct use of Wigner’s rational approximation. The substitution of Eq. 5.113 into 5.108, leads to the well-known “equivalent relation” whereby the flux exhibits the same functional form as that for homogeneous media, †.0/ p C †e ;

0 .u/ D .0/ †p C †e C †.0/ t

†e D

1 l

(5.117)

where †e is referred to as the “equivalent” cross section. It follows that .RI /k D .eq/

where p

p.eq/ a J.ˇk0 ; k ; ak / cos 2 l E0

(5.118)

i. h .0/ D †p C †e N0 is commonly referred to as the “equivalent” po-

tential scattering cross section and ˇk0 D p.eq/ =.0k cos 2 l /. Hence, Eq. 5.118 is identical to the corresponding resonance integral for homogeneous media defined by Eq. 5.81 if p.eq/ is replaced by p . It was the early version of this equivalence relation which was later generalized to the closely packed lattice as one shall see. This approximation was used extensively not only in reactor physics calculations but was also used as a tool to analyze resonance integral measurements. If one neglects the temperature dependence and the asymmetric component of the resonance line shapes, Eq. 5.118 for a low energy resonance immediately reduces to the following form, 1=2 a 0k 0 ı ˇ .1 C ˇ 0 / (5.119) .RI/k D 2 E0

262

R.N. Hwang

In the limit of a strong resonance, the above expression can be approximated by 1=2 S0 1=2 S0 .RI/k const †.0/ C / A C B p 4V0 M

(5.120)

which was used as the means to estimate the approximate behavior of the resonance integral as a function of the surface to mass ratio in analyzing experimental results. Similar results were also obtained by Gourevich and Pomeranchouk [7] based on Eq. 5.116. NRIM (or WR) Approximation In contrast to the NR-approximation, Eq. 5.108 under the NRIM-approximation in the same context described for homogeneous media must assume the following form: i h .0/ .0/ †t 0 .u/ D .1 P0 / †.0/ .u/ † s m 0 .u/ h i .0/ C †.0/ (5.121) CP0 †.0/ t .u/ †m m Again, like the case of the NR-approximation, one practical method to circumvent the task of computing the escape probability rigorously is the use of the rational approximation. Upon substituting Eq. 5.121 into 5.108, one obtains the pertinent equivalence relation as before, i.e., N †.0/ m C 1=l

0 .u/ D .0/ .0/ †m C 1=lN C †a .u/

(5.122)

Thus, the corresponding resonance integral for a Breit–Wigner resonance becomes .eq/

.RI/k D

p t J.ˇk0 ; k /; E0

.0/

p D

†m C †e N0

(5.123)

. where ˇk0 D p.eq/ .0k ak =tk / and the equivalent cross section is the same as that given in Eq. 5.117. It is quite clear that the same rationale is equally applicable to the IR approximation described in the previous section.

5.5.2 Traditional Collision Probability Treatment in a Closely Packed Lattice The discussions so far have been focused on the idealistic case in which the fuel lumps are considered as isolated from one another, separated by large moderator

5


263

regions. The basic principle was later extended to the more realistic situations involving closely packed reactor lattices. 5.5.2.1 General Features of the Escape Probability In the presence of many closely packed unit cells, neutrons escaping from the fuel lump in a given cell may not necessarily suffer their next collision in the surrounding moderator. This issue was first addressed by Dancoff and Ginsburg [60]. The “shadowing” effect attributed to neighboring fuel lumps on the escape probability of the fuel lump in question has henceforth been referred to as the Dancoff effect. The most comprehensive discussions on this subject are believed to be those given by Rothenstein [61], Nordheim [52], and Lukyanov [10]. The physical attributes of this phenomenon can be best illustrated by following a straight-line trajectory of the neutron path through the fuel lumps and the moderator regions. Let subscript 0 denote the fuel lump in question and 2; 4; : : :; 2N denote all neighboring fuel lumps, while 1; 3; 5; : : :; 2N C 1 denote the moderator regions sandwiched between the fuel lumps. If ln is the chord length of the path in the n-th region, and the local transmission probability between two surface points is denoted by n , then for the isolated fuel lumps previously described, one has .iso/

P0

D

1 †.0/ t l

.1 h0 i/

(5.124)

To account for the Dancoff effect, one must take into account the probability that neutrons survive without suffering any collision along the trajectory through the successive regions. Let P0 be the improved escape probability for the lump 0 in such a lattice. It is not difficult to see that P0 is expressible symbolically in the general form as follows: P0 D

1 .0/ †t l

h.1 0 / Œ.1 1 / C 1 2 .1 3 / C i

(5.125)

Similarly, the escape probability for the surrounding moderator, P1 , can be obtained by interchanging the indices in the equation above. However, the escape probability in its general form is obviously too complicated for practical applications when resonance structures in cross sections are involved. One commonly used simplified assumption is that all unit cells are taken to be identical, a configuration referred to as an infinite lattice with repeated cells. Under such an assumption, it follows that h2n i D h0 i D T0 ;

h2nC1 i D h1 i D T1 ;

n D 1; 2; 3; : : : ; N

(5.126)

The repeated cell assumption leads to four approximations most frequently used in applications. In particular, the first two of these provide the theoretical basis for extension of Wigner’s rational approximation to accommodate the Dancoff effect in closely packed lattices.

264

R.N. Hwang

“Black” Limit Approximation If the total cross section of the absorber is assumed to be large, all higher-order terms in the square bracket of Eq. 5.125 diminish. In addition, if the average of the product is replaced by the product of the averages, one obtains, P0

1 †.0/ t l

.1 h0 i/.1 h1 i/ D P0.iso/ .1 C /

(5.127)

where C is referred to as the Dancoff correction factor depending on the cross section and geometric configuration of the moderator alone. Thus, the rationalfunction-based approaches derived for isolated lumps are equally applicable to closely packed lattices if the Dancoff factor is known.

Nordheim’s Approximation The fact that the absorbers may not always be considered as “black” provides the need for further improvement. Since the linear terms in Eq. 5.125 are generally most dominant, all higher-order products in the ’s can be adequately approximated by the respective products of the averages as pointed out by Nordheim [52]. Thus, upon substitution of the averages and noting the geometric nature of the series for the repeated cells, Eq. 5.125 can be expressed in the closed form P0 D

.1 h0 i/.1 h1 i/ .1 C /P0.iso/ D .0/ .0/ 1 h0 i h1 i †t l 1 .1 l†t P0 /C 1

(5.128)

Infinite Slabs The case of closely packed repeated cells consisting of infinite slabs is probably one of the few cases where the escape probability can be represented analytically without simplifying assumptions. It is quite obvious from Eqs. 5.111, 5.125 and Table 5.1 that the escape and transmission probabilities for this type of configuration are generally expressible as a single integral. If the absorber and moderator with the thicknesses t0 and t1 , respectively, are either in a periodic (or reflective) arrangement, Eq. 5.125 can be reduced to a single integral of manageable form as pointed out by Corngold [62] and Rothenstein [61], namely, P0

D

1 2†.0/ t t0

1 X 12 E3 ..n C 1/ 0 C n 1 / C E3 .n 0 C .n C 1/ 1 / nD0

2E3 ..n C 1/. 0 C 1 //

(5.129)

5


265

.1/ where 0 D †.0/ t t0 and 1 D †t t1 are the optical thickness of the fuel and of the moderator, respectively. E3 .x/ is the exponential integral of order 3. It is quite obvious that the escape probability of an isolated fuel plate defined by .0/ Eq. 5.124 corresponds to n D 0 with transmission probability h0 i D 2E3 .†t t0 /. On the other hand, under the black limit assumption, the Dancoff factor in Eq. 5.127 is simply C D 2E3 .†.1/ t t1 /.

Repeated Cells with Cylindrical Configuration The infinite lattice in this case can be pictured as an ensemble of unit cells each consisting of a fuel rod and a cylindricized moderator region. In contrast to the infinite slabs, one general complication here is that the resulting integral generally involves a double integral according to the basic geometric properties of the cylinder given in Table 5.1. As discussed by Takahashi [63] and also by Rothenstein [61], it is possible to derive the same type of expression as that for infinite slabs except that it involves a series consisting of double integrals. In particular, the inner integrals are identifiable with the Bickley–Nayler function [64] of order 3. From the perspective of practical applications, especially fast reactor applications involving an extremely large number of resonances over a large energy span, the method derived on this basis is far more difficult to use for routine calculations than that for the infinite slabs. Other practical alternatives will be discussed later.

5.5.2.2 Fine-Tuning of the Rational Approximation for Routine Applications As discussed earlier, the dramatic simplification of the lattice physics calculations for the case of an isolated cell was made possible via Wigner’s rational approximation, leading to the equivalent relation whereby the lattice physics calculations can be treated in the same context as those for homogeneous media. The extension to closely packed lattices is obviously possible via either the black limit approximation or Nordheim’s approximation discussed in the previous subsection. By substituting Eq. 5.113 into 5.128, one obtains, P0

ı 1 lN .1 C / D .0/ D .0/ ı †t .u/ C ŒS0 =.4V0 / .1 C / †t .u/ C 1 lN .1 C / ŒS0 =.4V0 / .1 C /

(5.130)

It is interesting to note that, physically, the Dancoff effect is equivalent to the reduction of the surface to volume ratio. As long as the moderator cross sections are taken to be constant, or C is independent of energy, the equivalence relation is clearly applicable as well. From the NR-approximation consideration, Eq. 5.130 leads to an equivalent cross section of the following form,

266

R.N. Hwang

†e D

1 lN

a.1 C / 1 C .a 1/C

(5.131)

where a, referred to as the Bell–Levine factor [65, 66], provides more fine-tuning of the rational approximation. Typically, a is set to 1.33 and 1.08 for cylindrical and slab geometries respectively. The same argument can be extended to the NRIM-approximation defined by Eq. 5.123. The second issue is how to compute the Dancoff factor simply at run time. One obvious choice is to apply the same rational approximation proposed by Bell [65], .1/

B D 1 C D

l 1 †t

1 C l 1 †.1/ t

(5.132)

N where the moderator cross section, †.1/ t , is taken to be constant and l1 here is the average chord length for the moderator region. Hummel et al. [67] have found that Eq. 5.132 can be significantly enhanced by the following modification: D 1 C D B C B4 .1 B /

(5.133)

The approximation above was found to be adequate and has been incorporated into the MC2 -2 code [54] for routine applications.

5.5.3 Connections to Resonance Integral and Multigroup Cross-section Calculations For closely packed repeated cells, there are two types of resonance integral-based methods widely used in reactor applications, which are summarized below.

5.5.3.1 Rational Approximation and Approximate Flux Based Approaches As discussed in the previous section, the direct equivalence between the treatment of the resonance integral in heterogeneous media and that in homogeneous media can be established to various degrees of sophistication when applied in conjunction with traditional approximations to the slowing-down equation. Like the previously described case of isolated fuel lumps, one only needs to redefine the “equivalent” cross section †e according to Eq. 5.131. For the NR-approximation, it is readily amenable to the J -integral approach discussed in Section 5.4 from which the multigroup cross sections can be obtained via Eq. 5.92. The viability of this type of approach has been demonstrated by the MC2 -2 code. The same logic applies if one chooses to use the NRIM or IR approximation.

5


267

The equivalence relation also makes possible the use of the Bondarenko scheme widely used for routine reactor calculations. For the case in point, the precomputed self-shielding factors at a given temperature as a function of p.eq/ defined by Eq. 5.131 are generated for various preselected groups and stored prior to their deployment at run time.

5.5.3.2 Nordheim’s Method Like the case of homogeneous media, Nordheim’s method using a numerical method for solving the slowing-down equation and the subsequent calculation of the resonance integral described in Section 5.4.3.2 is also applicable here. The only difference is that the slowing-down equation is given by Eq. 5.108 instead of Eq. 5.78. Like Eq. 5.78, 5.108 can also be further simplified by assuming the NR-approximation for the diluents in the fuel lump in the same context as the former [52]. For fuel lumps with a single resonance absorber, Eq. 5.108 becomes 1 Zu .r/ 0 † .u / 1 0 s e .uu / .r/ F0 .u0 /du0 A F0 .u/ D .1 P0 / @ 1 ˛r †t .u0 / u"r .r/ C P0 †t .u/ C †m 0

(5.134)

where the indices r and m denote the resonance absorber and the diluent in the fuel lump respectively. The above equation would become identical to Eq. 5.78 if one sets the escape probability to zero. Thus, by using the escape probability defined by Eq. 5.128, the same numerical algorithm described in Section 5.4.3.2 for computing the resonance integral can be deployed. It is important to note that the original numerical scheme was developed primarily for computing a well-isolated Breit–Wigner resonance for thermal reactor applications. The question will arise when the mutual self-shielding effects resulting from the presence of other resonances in the vicinity of the one under consideration, a scenario quite common in fast reactor physics calculations. Therefore, more rigorous methods were subsequently developed in order to cope with the overlapping of resonances in the slowing-down process.

5.5.4 Rigorous Treatment of Resonance Absorption in a Unit Cell with Multiple Regions and Many Resonance Isotopes Motivated by needs in conjunction with fast reactor development, rigorous methods for treating resonance absorption were developed for various types of calculations where accuracy is required. The case in point can be viewed as a natural extension of the rigorous treatment of the slowing-down equation in homogeneous media if

268

R.N. Hwang

the same hfg approach is used as the means to discretize the lethargy domain. For a unit cell with multiple regions and many resonance isotopes, generally it would require the numerical solution of a system of N slowing equations analogous to the coupled equations for the two-region cell defined by Eqs. 5.104 and 5.105, i.e., Vi Fi .u/ D

N X

Pij .u/Vj Sj .u/

(5.135)

j D1

where Pij .u/ is the collision probability for neutrons originating in region j that will make their next collision in region i and Sj .u/ is the scattering source from region j as defined in Eq. 5.106. In principle, the numerical solution of these slowing-down equations can be carried out in the same context of the hyper-fine group approach described in Section 5.4.4.3. The analogy to Eq. 5.102 can be readily established via the use of matrix notation. For a given hfg k, the collision rates and scattering sources in various regions of the cell can be represented as vectors. In the same context as Eq. 5.102 for the homogeneous case, the collision rate can be expressed as .in/ E E D .P1 R/1 S E E .in/ E (5.136) CR D P S C R CR ; CR where R is a diagonal matrix with elements Ri i D .†sk =†t k /i .KQ s /i for a given region i and the remaining quantities appearing in Eq. 5.136 are consistent with Eq. 5.102 except that they are written in vector form. Thus, the collision rate for a given region can be specified once the collision matrix is known and the corresponding neutron flux is simply, CRi .uk / (5.137)

i .uk / D Œ†t .uk /i Given the hfg flux in every subregion of the unit cell, the desired broad group cellaveraged cross section for reaction process x required for practical applications becomes, Q xG

D

X N X i D1 k2G

x .uk / i .uk /

, X N X

i .uk / ;

x2i

(5.138)

i D1 k2G

The practical issue here is how to compute the collision probabilities for cells involving many regions and resonance isotopes at a large number of mesh points (or hfgs). For our purpose here, two proven methods will be outlined.

5.5.4.1 Kier’s Method for Cylindrical Unit Cells Consider a typical unit cell in the form of an infinite cylinder consisting of a fuel rod with cladding surrounded by a cylindricized moderator region [57]. Let us subdivide

5


269

the cell into N intervals all of which are in the form of annuli except for the interval at the center. To compute the detailed neutron flux at a given hfg and spatial interval requires an explicit or implicit description of the collision probability. This can be accomplished in steps as follows. Intra-Cell Neutronic Balance One way to simplify the complexity of computing the collision probabilities for the multiregion cell is to cast the intra-cell neutronic balance into a somewhat different perspective. In Kier’s approach [57], the interface current can be described in the form of a matrix equation, E TJE D PS (5.139) where the elements of T and JE denote the transmission probabilities and currents at the surfaces of each region, respectively, while P denotes the escape probabilE is the neutron source ities from a given region to its immediate neighbors and S in the region considered. Within the context of the hyper-fine group approach described in Section 5.4.4.3, this matrix equation must be solved at each hyper-fine group level with the source in each region computed the same way as before starting at the hfgs beyond the resonance region. Physically, one only needs to deal with the transmission and escape probabilities of the region under consideration to the immediately adjacent neighboring regions while the effects from the remaining regions are implicitly accounted for via boundary conditions. If one further assumes isotropic return of the currents at the boundaries, the corresponding T and P for annuli can be specified readily. To determine the interface currents based on the above equation requires the solution of a tri-diagonal matrix equation of order 2N 2N reflecting the tri-diagonal nature of the matrix equation given by Eq. 5.139 where two currents per region for such a region are required. Thus, Eq. 5.139 can be expressed in terms of 2N linear equations as specified in [54, 57]. As proposed by Kier [57], this matrix equation can be solved efficiently via the use of the usual Gauss elimination procedure followed by backward substitution once the pertinent transmission probabilities and escape probabilities are known. Computation of the Probabilities of a Given Annulus Consider the i-th annulus bounded by circles with radii ri 1 and ri respectively. For an annulus, three first-flight transmission probabilities and two escape probabilities are required to specify the neutronic balance defined by Eq. 5.139. The former can be denoted symbolically by TiOI ; TiIO and TiOO signifying the transmission of neutrons from inner-to-outer, outer-to-inner and outer-to-outer surfaces of region i , respectively. The latter can be denotedby PiC and Pi signifying the fraction

270

R.N. Hwang

of neutrons escaping the outer and inner surfaces of region i , respectively. These probabilities can be defined in manageable forms if the cosine current assumption is used, i.e., Z Z 2 x 4 4 ˛i TiOO D cos Ki3 .bi cos /d D p Ki3 .bi x/dx (5.140) sin1 ai 0 1 x2 q Z p z 4 1 OI 2 2 2 ; TiIO Dai TiOI (5.141) Ti D dxKi3 1 ai x ai 1 x 0 1 ai where ai D .ri 1 =ri /; ˛i D .1 ai2 /1=2 ; bi D 2†t i ri and z D bi .1 ai /=2. Here, Ki3 .x/ is the well-known Bickley–Nayler function [64] of order 3. The corresponding outward and inward escape probabilities for the i-th annulus are respectively, PiC D

1 1 1 TiOO TiIO ; Pi D 1TiOI †t i l †t i l

(5.142)

where lN D 2.ri2 ri21 /=ri 1 . Hence, all probabilities are specified once TiOO and TiOI are known. For practical calculations, these transmission probabilities can be either precomputed and stored in two-dimensional tables or alternatively, computed at run time via a low-order quadrature described in [68].

Determination of Collision Probabilities From neutron-conservation considerations, the collision rate for a given region can be specified symbolically as .CR/i D 1 TiOO Ji C 1 PiC Si ; i D 1 OO D 1 TiOI JiC TiIO Ji C 1 PiC Pi Si ; 1 C 1 Ti

i D 2; N (5.143)

where i D 1 denotes the inner-most region with the configuration of a circle. In the original development of Kier [57], the corresponding flux for region i and hfg k is simply, i .uk / D .CR/i =†t i .uk /. To put the results in the same context of Eq. 5.136, the traditional collision probabilities can be readily deduced from Eq. 5.143. Physically, it is not difficult to rationalize that the collision rate for region i is identifiable with Pij if one sets Si D 1 and Sj D 0 in the equation above, i.e., Pij D .CR/j ;

Si D 1; Sj D 0 if j ¤ i

(5.144)

Thus, the hfg flux for all subregions can be determined via Eq. 5.137, from which the cell-averaged broad group cross section defined by Eq. 5.138 can be specified. This approach has been incorporated into the RABBLE-code [57] as well as the MC2 -2 code [54].

5


271

5.5.4.2 Olson’s Method for Unit Cell with Many Plates During the period of fast reactor development, a great deal of emphasis was placed on the analysis of measurements from various fast critical assemblies made of drawers with a large number of plates. A much more rigorous alternative for such an approach was developed by Olson [58] without resorting to the assumptions of the cosine or cosine-square for the interface current made by Kier [57], conditions plausible for cylindrical cells. The rationale was that the neutron transport in the infinite lattice consisting of repeated cells with infinite slabs can be specified directly by the integral transport equation without resorting to the use of Eq. 5.139 as was evident from the earlier work of Corngold [62] and Rothenstein [61]. In the context of the hfg approach described previously, the scattering source for a given hfg k and i-th region, Si.k/ encompasses the sum of scatter-in and self-scattering terms described in Section 5.4.4.3. In the approach proposed by Olson [58], the source is allowed to vary linearly within each plate and is then used in conjunction with the integral transport approach to compute the intraand inter-cell currents. The resulting current at mean free path ; D †t i x, beyond the plate with optical thickness i can be expressed as E i / D J. ;

Si.k/ ŒE3 . / E3 . C i / C f . ; i / 2 i

(5.145)

where the first term is identifiable with the traditional result using constant source and the second term amounts to a correction term taking into account the linear variation of the source. f . ; i / involves a linear combination of exponential integrals of both orders 3 and 4 [58]. If a plate j is at mean free path away from the plate i , the collision rate in j due to the source in i is E E C j ; i / E i / J. CR.i ! j / D J. ;

(5.146)

For the infinite repeated cells with periodic plate arrangement, the contribution of all plates i to plate j become E 1 .i ! j / D CR

1 h X

E C j C mh; i / E C mh; i / J. J.

i (5.147)

mD0

where h is the optical thickness of the unit cell. A similar expression can also be derived for cells with reflective boundary condition. The infinite sum involving exponential integrals of order n 3 can be evaluated readily via the use of the Gauss quadrature developed for this purpose [58]. Thus, the collision probability Pij for a unit cell is expressible as h i E 1 .j ! i / E 1 .i ! j / C ŒCR Pij D CR Si

(5.148)

272

R.N. Hwang

By substituting Eq. 5.148 into 5.136 followed by inverting the resulting matrix, the flux in a given plate for the hfg k can be determined. The substitution of the flux so obtained into Eq. 5.138 yields the desired broad group cell-averaged cross section for applications. The approach described here was first coded in the RABIDcode [58] and later incorporated into the current MC2 -2 code [54] as a valuable option.

5.6 Treatment of Unresolved Resonances The treatment of resonance self-shielding effect in the unresolved energy range constitutes one of the most important topics in resonance theory especially in fast reactor applications. The range under consideration is usually defined as the energy span between the upper boundary of the resolved energy region and an upper bound beyond which the self-shielding effect becomes unimportant. For major actinides, it spans the range from the low keV region up to about 100 keV. This is the energy region where the Doppler-width becomes much larger than the resolution width of the instrument so that resonances may not be parameterized deterministically. The methods for treating resonance absorption in this range can be viewed as a natural extension of the statistical theories of average cross sections such as those described by Moldauer [69] and Ericson [70]. Since the theoretical foundations may often be obscured in routine applications, it is useful to summarize briefly some conceptual aspects of the problem prior to the discussions of the basis for the computational methods. For our purpose here, the discussions will be focused on current methods based on the single-level Breit–Wigner approximation although the same concept is extendable to the more rigorous S -matrix approach if its parameters are expressed in terms of the R-matrix parameters.

5.6.1 Statistical Theory Basis The statistical basis for evaluating the average cross section hx iE is assumed to be extendable to the estimation of the expectation values of the reaction rate hx iE;Er and the flux h iE;Er as well. Without loss of generality, each microscopic cross section is represented by the R-matrix formalism defined previously in terms of parameters E and c . From the statistical theory of spectra [71], the distributions of these parameters are well-known. These distributions, in effect, define the joint density function (p.d.f.) required for evaluating the averages attributed to an ensemble of resonances ˝ 2in˛ the vicinity of an energy, say E . Given information of hjEi Ei C1 ji and ci and through the explicit knowledge of the behavior of x and as function of energy, the expectation values of interest are, in principle, completely specified.

5


273

5.6.1.1 Some Statistical Theory Fundamentals A brief outline of some pertinent statistical theory fundamentals is believed to be a good starting point for the case in point.

Basic Rule The probability distribution of an event A characterized by statistically independent variables .x1 ; x2 ; x3 ; : : : ; xi / 2 A is defined as Z Z P .A/ D

:::

Z Y

pi .xi /dx1 dx2 : : : dxi

(5.149)

i

where pi .xi / is referred to as the probability density function (p.d.f.) for xi and the products of all p.d.f. is known as the joint density function for these variables. For the case under consideration here, these quantities are identifiable with the partial width distributions and the level correlation functions. The direct use of the joint density function provides the most widely used basis for computing the averages of interest to be described.

Addition Theorem The probability distribution function for the union of two events is given by, P .A [ B/ D P .A/ C P .B/ P .AB/

(5.150)

Multiplication Theorem The probability distribution of the product of two correlated events is given by, P .AB/ D P .A/P .BjA/ D P .B/P .AjB/

(5.151)

where P .BjA/ is known as the conditional probability of B for a given occurrence of the event A. The above relation must be symmetric when A and B are interchanged. This theorem is of a great deal of practical interest since it provides another basis for computing the averages of interest known as the “probability table” method to be described shortly.

274

R.N. Hwang

5.6.1.2 Statistical Distributions of Practical Interest There are three types of distributions for the resonance parameters of practical interest. They will be summarized as follows.

Partial Width Distributions It is well-known that the reduced width amplitude c for a given level and channel c follows the normal distribution with zero means and variance of unity according to Porter and Thomas [72]. For a given reaction process, there may exist several channels and, consequently, the evaluation of the average quantities can require the task of evaluating many multiple integrals with many variables. For the Breit–Wigner approximation, the partial width of a given reaction type is actually used. For a given reaction process x, it is defined as X 2 2Pc c (5.152) x D c2x

where Pc is the penetration factor described in Section 5.2.1.1. If one assumes that 2 > is the same for all c 2 x, the probability density function for the average < c the partial width becomes, p .y/dy D

2 1 y e 2 y dy 2 .=2/ 2

(5.153)

where p .y/ is the well-known 2 -distribution with the degrees of freedom equal to the total number of exit channels for process x and y D x =hx i, the local to average ˝ 2 ˛ ratio of the partial width. It should be noted that Eq. 5.153 is no longer valid is different for c 2 x, a scenario where the multichannel effect is important. if c Level Width Distribution The distribution of E is characterized by the Wigner distribution [73] or the longrange correlation described by Dyson [74]. The former, also known as Wigner’s surmise developed in 1957, that the probability distribution of the level spacing jE E j for a given l and J state should follow a simple analytical form, w.x/dx D x exp x 2 dx (5.154) 2 4 ı˝ ˛ where x D jE E j jE E j is the local to average ratio of the level spacing. The fact that w.0/ D 0 signifies the repulsion of eigenvalues of a real and symmetric Hamiltonian matrix. More studies of this subject of such a matrix ensemble were further pursued by Wigner and others [71, 73, 74]. Of particular interest for

5


275

practical applications is the so-called Gauss orthogonal ensemble (GOE) in which the distributions of the matrix elements are taken to be statistically independent Gaussians and invariant under orthogonal transformation. The Wigner distribution given above is identifiable with the behavior of the difference of eigenvalues of a 2 2 matrix. In contrast, the latter proposed by Dyson [74] is free from the assumption used by the former. Here, the ensemble under consideration is not the Hamiltonian directly but an ensemble of eigenvalues fE g of a unitary matrix S of the form exp.i ‚ / uniformly distributed around a unit circle. For an orthogonal ensemble, a general expression for the n-level correlation function as a function of ‚ with D 1; 2; 3; : : : ; n was developed by Dyson [74]. From the practical point of view, the general form of the distribution is difficult to use. In the limit of two levels, however, it is expressible in a simple analytical form which can be utilized for routine applications as one shall see. Level Correlation Functions In practical calculations, another distribution of interest, .jEk Ej j/, is the probability of finding any level within an interval at a distance jEk Ej j away from a given level Ek when the direct integration approach is used for averaging. It is of particular interest when the overlapping ˇ˛ is dominated by the immediate ı˝ˇ of levels neighbors. If one sets x D jEk Ej j ˇEk Ej ˇ as before, this distribution can be defined by the following convolution integral equation, Zx .x/ D w.x/ C

w.t/.x t/dt

(5.155)

0

where w.x/ is the Wigner distribution of the level spacing belonging to a resonance sequence of a given l and J previously defined. The analytical solution to this equation in a closed form is difficult to derive but one can handle it numerically. It is interesting to note that the analytical solution in a closed form can be derived if w.x/ is replaced by the 2 -distribution of even order via the use of the Laplace transform [75], i.e., 1 .x/ D 2 i

cCi Z 1

ci 1

.v=2/v=2 e xp dp .v=2/ .v=2/x e D 2 i .p C v=2/v=2 .v=2/v=2

cCi Z 1

ci 1

e .v=2/xz dz z.v=2/ 1 (5.156)

whereby, for all even 2, it is expressible in terms of a damped oscillatory function in closed form readily shown via the use of the Cauchy integral theorem with the poles of the integrand computed via De Moivre’s theorem. In particular, by setting D 8, the corresponding 2 -distribution exhibits resemblance to the Wigner distribution and Eq. 5.156 yields the following simple analytical form, .x/ D 1 e 8x 2 sin 4xe 4x

(5.157)

276

R.N. Hwang

which was used in earlier work based on the “nearest level” approximation for estimating the overlap effects on the self-shielded cross sections attributed to immediate neighbors. Physically, Dyson’s two-level correlation function can be considered as the rigorous alternative to .x/. The former was later introduced into reactor applications by replacing .x/ with r2 .x/ derived by Dyson [74], i.e., ds.x/ r2 .x/ D 1 Œs.x/ dx

Z1

2

s.y/dy; s.x/ D

sin x x

(5.158)

x

There is another type of level correlation function of practical interest. If the ensembles levels fEk g and fEj g belong to either different spin states or different nuclides, the distribution of jEk Ej j must follow the Poisson distribution. Thus, if one substitutes the 2 -distribution of D 2 into Eq. 5.156 in place of the Wigner distribution, the corresponding correlation function becomes .x/ D 1. Physically, this signifies the fact that the probability of finding a level at a distance jEk Ej j is equally probable so long as these levels are statistically independent.

5.6.1.3 Conceptual Aspects of Computing Average Cross Sections Conceptually, the computation of these averages can be considered as the natural extension of the previously described treatment of resonance cross sections in the resolved energy range. With no loss of generality, consider an ensemble of a generic quantity qk.l;J / .E/ attributed to a Breit–Wigner resonance k for a given l and J state. E D , can be pictured as the population The statistically averaged value, qk.l;J / .E/ E par:

average in the following context, D

E

qk.l;J / .E/ E par:

N Z E2 X 1 D qk.l;J / .E/dE E2 E1 E1 kD1 Z 1 1 1 .l;J / lim D q .E/dE hDi N !1 N 1 k Z 1 1 q .l;J / .E/dE D hDi 1 k par:

(5.159)

where the average resonance parameters for the distributions are specified at the predetermined energy E provided by the evaluated data file. Here, hDi denotes the average level spacing. The equivalence of the discrete average and the statistical average clearly requires the validity of the ergodicity and stationarity of the samples inside the angular bracket in the vicinity of E .

5


277

This concept has been widely used as a basis for computing the unshielded, as well as the shielded, group cross sections. For the former, qk.l;J / .E/ is set to .l;J / .E/ as defined in Eqs. 5.41 and 5.42. For the latter, it requires the average rexk .l;J / action rate as well as the average flux, which can be obtained by setting qk .E/ .l;J / equal to xk .E/ k .E/ and k .E/ respectively. The relation defined in Eq. 5.159 can be cast into a different context from the standpoint of the widely used Monte Carlo technique. Given the probability density functions for level spacing and partial widths, one can construct the corresponding discrete resonance sequence whereby the averages can be treated in the same manner as the resolved resonances as shown in the first part of Eq. 5.159. This can be accomplished by the following procedure. First, one specifies the c.d.f. from the known p.d.f., say p.x/, Zx P .x/ D

p.y/dy; x D P 1 .x/

(5.160)

0

where P 1 .x/ symbolically R 1 represents the value x obtainable from a given value in P .x/. By definition, 0 p.x/dx D 1 so that the range of the c.d.f. must be between 0 and 1. P .x/ can be either the c.d.f. for level spacing or that for partial widths where appropriate. Second, a resonance sequence can be generated if one assumes that P .x/ is equally probable within the range 0 P .x/ 1. This can be accomplished by choosing a random number i obtained by a random number generator, a computer program for generating statistically uncorrelated numbers with 0 i 1, and setting Pi D i , from which the corresponding set of parameters xi can be obtained from the inverse relation given by Eq. 5.160. Thus, by repeating such a process for the c.d.f.’s of all resonance parameters, a resonance sequence consisting of discrete resonances can be generated. A sequence of discrete resonances so generated is commonly referred to as a “resonance ladder.” Given the resonance ladder, the desired average can be treated in the same manner as if they were resolved resonances. It is important to realize, however, that this procedure is always accompanied by significant statistical uncertainties reflecting the highly fluctuating nature of the resonance phenomena. Thus, to reduce such uncertainties generally requires the inclusion of many statistically uncorrelated ladders. This undoubtedly hinders the routine usage of such a method especially when high accuracy is required. It should be noted that the desired averages can also be specified from a somewhat different perspective via the use of the multiplication theorem defined by Eq. 5.159. This concept leads to another alternative approach commonly referred to as the “probability table” method, which is particularly attractive in conjunction with reactor calculations using the Monte Carlo method. A discussion of this method will be presented later.

278

R.N. Hwang

5.6.2 Average Unshielded Cross Sections and Fluctuation Integrals The simplest possible averages based on the statistical concept given in the previous subsection are the average unshielded cross sections also referred to as the “infinitely dilute” cross sections. Such averages are generally temperature-independent due to the invariant nature of the Doppler-broadened line-shape functions when integrated over energy. One average of particular importance is the average compound nucleus cross section. Substituting Eq. 5.42 into 5.159, one obtains, htR i D 2 2 2

X

gJ cos 2 l sl;J

(5.161)

l;J

where sl;J D hn i=hDi is known as the “strength function,” a quantity directly relatable to the transmission measurements. For partial cross sections, the unshielded averages are directly expressible in terms of fluctuation integrals. The substitution of Eq. 5.41 into 5.159 leads to, X g n x J ; x 2 ; or f hx i D 2 hDiJ t l;J X g 2 n J hsR i D 2 2 2 2 sin2 l hn i hDiJ t 2 2

(5.162)

(5.163)

l;J

for reaction x cross section and resonance scattering cross section, respectively. It should be noted that the total width in practical applications is taken to be the sum of all partial widths plus a “competitive width” c , i.e., t D n C C f C c . The latter presumably provides the approximate means to account for the inelastic scattering channels in the energy range above the first excited state of the nuclide in question. The quantity hn x =t i is referred to as the “fluctuation integral” which can be computed accurately via direct integration over all distributions of the partial widths. Let fx D hn x =t i be the multiple integral over the 2 -distributions of partial widths involved. Taking advantage of the exponential nature of these distributions, one can show readily that fx can be reduced to a single integral of the generic form given below, Z1 fx D aC 0

e t dt .1 C c1 t/. 1 =2/Cj .1 C c2 t/. 2 =2/Ck .1 C c3 t/. 3 =2/Cm

(5.164)

for x 2 ; s, or f , where ci D 2hx i = = i . with the degrees of freedom ˝i in ˛ the context given by Eq. 5.153 and i 2 1; 2; 3 corresponding to hx i D hn i ; f , and hc i, respectively, and C D hn i hx i = . The constant a is equal to unity for all x ¤ s while a D 3 when the fluctuation integral for the scattering cross section

5


279

that appears in the first term of Eq. 5.163 is considered. The integer constants j , k and m are dependent on the type of reaction considered. For x 2 , j D 1, k D 0, and m D 0. For x 2 f , j D 1, k D 1 and m D 0. For x 2 s, j D 2, k D 0 and m D 0. For x 2 c, j D 1, k D 0 and m D 1. It should be noted that fc is not used in practical calculations. One way to evaluate the single integral given by Eq. 5.164 is via a Gauss quadrature. In particular, the integral can be computed accurately and efficiently if the change of variable of integration setting y D .1 C t/=.1 t / is made prior to integration and is followed by the use of the Gauss–Legendre quadrature in the y-domain [77]. The accurately computed average unshielded cross sections via the fluctuation integrals provides the necessary criteria to verify those computed by other methods, either via deterministic or Monte Carlo approaches. The disagreement in the unshielded cross sections clearly will be passed on to the calculations of the selfshielded cross sections as well.

5.6.3 Traditional Approaches for Computing Self-shielded Average Cross Sections Discussions in this section focused on the traditional methods for computing the self-shielded cross sections in the unresolved resonance range that directly utilize the known joint probability distribution of resonance parameters described in Section 5.6.1.1. There are two different types of approaches based on the same general principle currently in use. The most commonly used approach is based on the direct integration, a procedure similar to that of computing the unshielded averages. Another less commonly used approach, referred to as the “resonance ladder” method, is based on statistically generated discrete resonance sequences via the use of the Monte Carlo technique described in Section 5.6.1.3. A brief discussion of these methods is now presented.

5.6.3.1 Methods Based on Direct Integrations For the direct integration approach, the self-shielded cross section for a given reaction process x associated with a given l and J state is defined on a statistical basis in the vicinity of a given energy E and can be viewed as the natural extension of the J -integral approach described in Section 5.4.4.2 if the NR-approximation is assumed [53]. The only difference is that the sum over individual resonances is replaced by a statistical average in the context of Eq. 5.159, i.e., D Q x.l;J / D

†b

x.l;J / .E / †b C†Rt .E /

f

E

; f D1

†Rt .E/ †b C †Rt .E/

(5.165)

280

R.N. Hwang

where †b denotes the macroscopic energy-independent background cross section including all potential scattering cross sections and smooth cross sections given by the data file as well as the “equivalent” cross section for the cell where appropriate. Like the case of the resolved resonances, the presence of a neighboring resonance can affect the resonance integral of a given resonance through its impact on the local flux. From statistical considerations, such overlap effects can come from two types of resonances present. They may either belong to resonances of the same spin sequence or those of a different spin sequence of the same nuclide and/or sequences of different nuclides. The former requires statistical averaging over the level correlation function defined by Eq. 5.157 or 5.158, while the latter requires integration over the uniform distribution, or .x/ D 1, as described in Section 5.6.1.2 in addition to averaging over the partial widths. Two scenarios can be described separately as follows.

Scenario 1: One Spin Sequence Only In absence of other statistically independent resonance sequences, Eq. 5.165 can be expressed as [53],

Q x.l;J /

h iE x J.k ; ˇk ; ak / Ok.x/ 0k D h iE D 1 1 hDi t J.k ; ˇk ; ak ; ak / Ok.t0 /k b hDi

D

(5.166)

where b D †b =N0 is the background cross section per atom density of the nuclide under consideration. As demonstrated in [53,77], the statistical average can be computed efficiently via the use of Gauss quadratures. In particular, the integration involving the level correlation function can be significantly simplified if the overlap effects can be attributed to the nearest level, an approximation usually referred to as the “nearest neighbor” approximation.

Scenario 2: Presence of Many Statistically Independent Resonance Sequences As described in [53], the integrand in Eq. 5.165, in principle, can also be further partial-fractioned into the similar form given by Eq. 5.166 in the presence of other statistically independent resonance sequences. Let the subscript i denote such a sequence and Jk denote the corresponding J -integral. The statistical averaging of such an integral requires not only integration over the correlation function of neighbors of the k-th sequence but also a uniform distribution of those belonging to the i-th sequence as well. The latter integration to first order can be approximated readily and the resulting expression for the reaction rate and the flux become,

5


281

2

E E X 1 b D b D .x/ .x/ 4 xk Jk x Jk 1 hDk i hDk i hDi i

3 E D .t / 5 ti Ji (5.167)

i ¤k

f 1

X all m

E 1 D tm Jm.t / hDm i

(5.168)

where Jk.x/ denotes the quantity inside the square bracket in the numerator of .t / Eq. 5.166 while Jk denotes the quantity inside the square bracket in the denominator of Eq. 5.166. Note that the second-order terms in Eqs. 5.167 and 5.168 are usually small and the inter-sequence overlap terms can be approximated by, E Y E X 1 D 1 D .t / .t / 1 ti Ji ti Ji 1 hDi i hDi i all i

(5.169)

all i

Thus, the substitution of Eq. 5.169 into 5.167 and 5.168 leads to an extremely simple expression for the average self-shielded cross section particularly useful for reactor applications, X Q x.l;J / (5.170) Q x D l;J

where Q x.l;J / is dependent on the in-sequence overlap term but independent of the inter-sequence overlap terms as given by Eq. 5.166. Physically, the direct integration approach is directly compatible with the fluctuation integral approach for computing the unshielded average cross sections. In the limit of infinite dilution, the average shielded cross section will approach the unshielded cross section based on the fluctuation integral.

5.6.3.2 Methods Using the Statistically Generated Resonance Ladders As described in Section 5.6.1.3, the unresolved resonances can also be pictured as an ensemble of discrete resonances generated based on the known distribution functions of resonance parameters via the usual Monte Carlo techniques. One advantage of this approach is that subsequent calculations of the shielded average cross sections are no longer constrained by the use of the NR-approximation required by the previous methods based on direct integration. However, the resulting cross sections computed on this basis are subject to large statistical uncertainties difficult to circumvent. Therefore, this approach is not widely used in practical applications as others.

282

R.N. Hwang

5.6.4 Probability Table Methods One method particularly amenable in conjunction with the Monte Carlo approach for reactor applications is the “probability table” method pioneered by Levitt [78]. The method bypasses the need of directly using the ladder approach at run time in earlier Monte Carlo codes, a procedure that would not only require large storage space and excessive computing time but also would lead to less manageable statistical uncertainties in its results. Conceptually, the general idea is to utilize the probability and conditional probability distributions themselves deducible numerically from the known distributions of resonance parameters so that the averages of interest can be cast into the form of a Lebesgue integral instead of that of the Riemann integral described earlier.

5.6.4.1 Conceptual Basis Consider a simple case where the neutron flux .t / is a function of t alone, i.e., the NR-approximation. One alternative joint probability distribution required to specify the average of the type hx .t /i is, according to the multiplication theorem, .max/ Z t

hx .t /i D

.max/ Z t

p.t /E.x jt / .t /dt ; h .t /i D 0

.t /p.t /dt 0

(5.171) where p.t / is the p.d.f. for the total cross section and E.x jt / is referred to as the conditional means of the partial cross section, the average x over the conditional probability p.x jt /, defined as .max/ Z x

E.x jt / D

x p.x jt /dx

(5.172)

0

which can be precomputed once the conditional probability is known. Thus, if these probability distributions are deducible from the known distributions of resonance parameters, the averages of practical interest can be expressed in terms of single integrals of Lebesgue form given by Eq. 5.171. For practical applications, only the prior knowledge of p.t / and the conditional means of various partial cross sections is required at run time. The main attraction of this approach is that p.t / along with its c.d.f. and various conditional means can be predetermined and tabulated in one-dimensional arrays for a prescribed E and different temperatures before their deployment. The same principle is extendable to the treatment of the resolved resonances in conjunction with the subgroup methods developed

5


283

independently by Nikolaev [79], Cullen [80] and Ribon [75]. Various numerical methods for computing these tabulated quantities will be further addressed. There is one word of caution if this approach is applied beyond the NR-approximation. Strictly speaking, this would require the conditional distribution of the form p.x js ; t /, which may negate, at least in part, the advantage of simplicity described above.

5.6.4.2 Methods for Computing the Tabulated Quantities There are two widely used methods for computing the discrete values of p.t / and E.x jt / as a function of t . One is the numerical scheme based on the statistically generated resonance “ladders” originally developed by Levitt [78]. The other is the moment-based approach based on matching the moments and partial moments of cross sections developed by Ribon and Maillard [75]. A brief outline of these methods is presented below.

Method Based on Resonance Ladders The rationale is relatively straightforward. The discretized p.t / and the corresponding conditional means are determined through the following steps. 1. For a given reference energy E where the averages are to be computed, a set .max/ is preselected as the of total cross section bands fti g with 0 < ti < t mesh points for the tables to be constructed. t D ti1 and t D ti , graphically equivalent to two parallel lines in a plot of t vs. E, constitutes the i-th “band.” The index i will henceforth be used to denote the t bands that will serve as entries to the tables. 2. Select an energy range E D E2 E1 around E with E >>> hDi, the average spacing of the resonance ladder in question. 3. Generate a “ladder” of resonances with the peak resonance energies falling inside E using the usual Monte Carlo scheme described in Section 5.6.1.3. 4. Construct a set of fine energy mesh fEk g within the preselected interval E from which the corresponding point-wise total cross sections ftk g and partial cross sections fxk g are computed. 5. Determine a subset of energy pairs fE2j ; E2j 1 g, signifying the consecutive pairs of intersections created between the upper and lower lines of the i-th “band” with the plot of t vs. E, which can be obtained via the criterion ti < tk < ti1 . The index j will henceforth be used to denote the consecutive energy pairs sandwiched by a given band i . 6. Store the localized average partial cross sections N xj.L/ .E2j ; E2j 1 / and total .L/

cross section N tj .E2j ; E2j 1 / of this ladder L for all energy pairs within the given band.

284

R.N. Hwang

7. The localized average of the probability density function (p.d.f.) of the total cross section at the i-th band amounts to the fraction of N tj.L/ .E2j ; E2j 1 / that falls in the range between ti and ti1 , i.e., J P

pi.L/

D

j D1

jE2j E2j 1 j E

D 0;

; ti1 < .L/ tj .E2j ; E2j 1 / < ti elsewhere

(5.173)

if and only if the band width is sufficiently small. .L/

8. The c.d.f. for the given ladder is simply Pi

D

i P

kD1

.L/

pk .

9. The corresponding conditional means of the partial cross section at the i-th band becomes E.x j ti / D .L/

1 X .L/ .L/ .E2j E2j 1 / xj .E2j ; E2j 1 /; ti1 .v†/min that they observed, specifically in small heavy water assemblies and small beryllium assemblies. This apparent paradox was subsequently resolved by analysis using simple cross-section and scattering kernel models – that retained the salient physical features – in conjunction with the exp.i BE rE/ spatial ansatz. Because of the simple models used, the analysis – done using a simple Laplace transform in time – could be carried out rather explicitly, and it showed that for a detector response the continuous spectrum or singular region, shown shaded in Fig. 8.4 leads to a branch cut (see Fig. 8.5), and that, for a sufficiently small system (large B 2 ), where the fundamental eigenvalue ˛0 (pole in the Laplace transformed solution) has passed through the branch point at .v†/min , it bifurcates, and the two poles that are born in this bifurcation move onto the adjacent sheets of the Riemann surface (see Fig. 8.5). The resulting time dependence for the detector response, given by the integral around the branch cut, is very close to exponential for a long time (but not as t goes to infinity) if the bifurcated poles are very close to the branch cut [48]. Of course, if this pseudo-exponential decay is very close to exponential and if it lasts for a long time, it would be observed in a pulsed experiment as exponential decay – since, by the time (slower) non-exponential behavior would appear the count rate would most likely be so low as to be indistinguishable from the background count in the experiment. A separate, not unrelated analysis, based on a more specialized cross-section model for a crystalline material with a Bragg peak (e.g., beryllium), gave an earlier somewhat different explanation of exponential decay rates above .v†/min in crystalline materials [49]. ® ® s = −iv . B − vΣT (v)

Im (s)

s = −min[vΣT (v)]

· · · - a2

- a1

- a0

Re (s)

Fig. 8.4 Sketch of the complex s-plane for the velocity-dependent transport operator in the image reactor theory approximation

410

J. Dorning

Im (s)

Branch Cut X X

(nS)min

Re (s)

Bifurcated Poles on Adjacent Sheets of the Riemann Surface

Fig. 8.5 Sketch of the complex s-plane for the energy-dependent diffusion operator

Many important theoretical and experimental problems in neutron thermalization related to thermal reactor physics – not just those related to pulsed neutron experiments – were put to rest during those years. Cross-section measurements and time-of-flight measurements of energy and angular distributions in thermal neutron scattering experiments, combined with the development of theoretical models for slow-neutron scattering, led to very much improved cross-section and scattering kernel data for use in reactor design and analysis. (See [24] for a very nice summary of some of this research written shortly after the research was completed.) These data were the antecedents (and in many cases the origins) of the Evaluated Nuclear Data Files (ENDF), and many of them are still used today. These were exciting times! They were great fun! And some really good science and engineering was done along the way! And the present author was thrilled to be part of it – albeit only at the very end. Like great wines, some data improves with age! (And even becomes perfect!) Before leaving the subject of pulsed neutron experiments in thermal systems – now so distant in time, but still so close to my heart – a personal anecdote seems irresistible (and perhaps even justified, since it was part of the original oral presentation of this material). Figure 8.6 from [50] shows experimental values (all well below .v†/min / of the decay constant in pulsed neutron experiments in light water assemblies along with dashed and solid curves that represent the calculated values of ˛0 based on the Nelkin model [51] and its slight improvement, the so-called anisotropic model [52]. It is clear that the data point for the smallest experimental system – actually a pingpong ball into which water was injected using a hypodermic needle – does not lie on the solid curve. And in the original graph given to the draftsman (by me), it was even oh so slightly further from the curve. (It took a few iterations to get where it is on this figure. After that I gave up!) Thus, Fig. 8.6 shows the initial improvement of that data point with age. But the real improvement was yet to come – when, soon

8

Nuclear Reactor Kinetics: 1934–1999 and Beyond

411

10 9 8

α0 (sec−1), ×104

7 6 NELKIN MODEL 5 4 3 ANISOTROPIC MODEL

2 1 0

0

1

2

3

4

5

6

7

8

9

10

11

R (cm)

Fig. 8.6 Experimental points and theoretical curves for the decay constant in spheres of water as a function of their radius (with permission from the American Nuclear Society, Copyright July 1968 by the American Nuclear Society, La Grange, IL)

after [50] appeared, the textbook by G. Bell and S. Glasstone [7] was published with an apparently redrafted version of this figure in it, on which this data point lies smack dead center on the solid curve! “Like great wines. . . .!”

8.4.3 Exponential and Non-exponential Decay in Subcritical Fast Multiplying Assemblies The simplest description of a pulsed neutron experiment performed in a subcritical, multiplying assembly rather than a nonmultiplying system is again the time-dependent one-speed diffusion equation, Eq. 8.65, but with v†a replaced by .1 ˇ/k1 v†a . Then the same simple steps reviewed at the beginning of Section 8.4.1 lead to ˛0 D v†a Œ1 .1 ˇ/k1 C vDB 2 ;

(8.85)

from which it immediately follows that keff

k1 D 1 C L2 B 2

1 ˛0

`1 1CL2 B 2

.1 ˇ/

D

1 ˛0 `0 ; .1 ˇ/

(8.86)

412

J. Dorning

where, as indicated below Eq. 8.6, `1 D 1=v†a , etc. and the overbar has been added here as a reminder that the v†a which appears in the one-speed equation has been averaged over the neutron energy spectrum. It is clear from Eq. 8.86 that a pulsed neutron experiment in a subcritical multiplying assembly is an integral experiment that yields the effective multiplication constant keff for the assembly; and from Eq. 8.85 that a sequence of such experiments on successively larger assemblies (smaller values of B 2 ) gives k1 for the system, given v†a and ˇ are known. Not unexpectedly, when the velocity-dependent Boltzmann equation is used to describe even a subcritical thermal multiplying assembly, all the subtleties summarized in Section 8.4.2 above arise. And in a subcritical fast multiplying assembly they are more pronounced. This is because, unlike in a thermal assembly in which the up-scattering of the thermal neutrons provides an energy-spectrum regeneration mechanism (in addition to that which results from the fission process) that tends to lead to the establishment of a collective energy mode, in a fast assembly, in which only down-scattering occurs, the resulting downward shift of the energy spectrum that occurs continuously in time does not lead to such a mode. Only when a fast assembly is quite close to critical, and the fission regeneration is dominant over the slowing down, is a persistent (decaying) collective energy mode established. From the mid-1960s to the mid-1970s, when various Western European countries, the Soviet Union, Japan, and the United States were developing plans to build fast reactors, pulsed neutron experiments were done as integral experiments to measure keff and other integral parameters, and to verify computational capabilities and models. Because they required considerably more financial and other investment (subcritical fast multiplying assemblies with plutonium or enriched uranium fuel) than non-multiplying thermal assemblies (e.g., spheres of water or graphite stacks) a limited number of these pulsed neutron experiments were done. Experimental data from two such experiments [53], done at the SUAK facility in the mid-1960s at the Kernforschungszentrum in Karlsruhe, Germany, are shown in Fig. 8.7 (adapted from [54]). The die-away data for SUAK-B – the larger, closer to critical assembly – clearly is exponential; whereas that for SUAK-A – the smaller farther subcritical assembly – clearly is not. Both sets of data initially appear to be exponential (upper portions of the straight lines; however, while the data for SUAK-B remain on the straight line at later times, those for SUAK-A depart from the straight line indicating slower than exponential decay at long times. This corresponds to pseudo-mode decay in which the die-away appears to be exponential for a finite time but then becomes slower at later times when the pseudo-mode ceases to be sustained and collapses. Because the .v†/min limit and measurements of related non-exponential die-away in thermal non-multiplying assemblies were well known by the time these experiments in fast assemblies were carried out, the nonexponential die-away in SUAK-A was not a surprise – and some experimentalists were, indeed, very well prepared for it. Similarly, some theorists also were prepared to do the analysis. Because of the fast, rather than thermal energy spectra in these assemblies – and especially because of the continuously downward-shifting fast spectrum in the pseudo-mode die-away – a one-speed diffusion theory description certainly was inadequate.

8


413

A : 25 ELEMENTS, keff =0.78 EXPERIMENTAL POINTS I/α = 122 nsec

5

10

B : 36 ELEMENTS, keff =0.869 EXPERIMENTAL POINTS AND CURVE, I/α = 230 nsec

4

10

A

B

3

10

10

2

0.2

0

0.4

0.6

0.8

1.0

1.2

1.4

1.6

1.8

2.0

t = ( m sec)

Fig. 8.7 Experimental die-away curves for the two SUAK assemblies (Adapted from [54])

A speed-dependent (or energy-dependent) representation was essential. And, since these fast assemblies were larger (even in mean free paths) than the most interesting thermal assemblies in previous studies, a speed-dependent diffusion theory description seemed adequate. Whether this description, or a speed-dependent transport theory description, is used it is simple and convenient to exploit the separability of the kernel in the Fredholm integral associated with the fission process. Starting from the speed-dependent diffusion equation with isotropic scattering (slowing down) in the center-of-mass system, introducing the ansatz ˆ.r; v; t / D '.B; v; t / exp.i BE rE/ and Laplace transforming from time to the variable s leads immediately to

s C vD.v/B C v†T .v/ 'Q B 2 ; v; s 2

Z1 d v0 v0 †s v0 ! v 'Q B 2 ; v0 ; s v

Z1 D .1 ˇ/.v/ d v0 v0 †f v0 'Q B 2 ; v0 ; s C 'Q B 2 ; v; 0 S B 2 ; v; s ; (8.87) 0

where '.B Q 2 ; v; s/ depends only on B 2 not B, and '.B Q 2 ; v; 0/ is the initial condition. Now, introducing the Green’s function of the associated slowing-down equation s C vD.v/B 2 C v†T .v/ GQ vjv0 ; B 2 ; s Z1 d v0 v0 †s v0 ! v GQ v0 jv0 ; B 2 ; s D ı.v v0 /; (8.88) 0

414

J. Dorning

which is simply the textbook slowing-down equation with v†t .v/ augmented by s C vD.v/B 2 , the solution to Eq. 8.87 then can be written in terms of this Green’s function as Z1

'Q B 2 ; v; s D

d v0 GQ vjv0 ; B 2 ; s S B 2 ; v0 ; s

0

Z1 D .1 ˇ/ d v0 GQ vjv0 ; B 2 ; s .v0 / Q B 2 ; s 0

Z1 C

d v0 GQ vjv0 ; B 2 ; s ' .B; v0 ; 0/ ;

(8.89)

0

where the Laplace transformed fission neutron production rate

2

Z1

Q B ; s

d vv†f .v/'Q B 2 ; v; s ;

(8.90)

0

has been introduced. Then multiplying Eq. 8.89 by v†f .v/ and integrating with respect to v from 0 to 1 yields a formal solution for the transformed fission neutron production rate, or the transform of the time dependence of a fission detector. Q B 2; s Q.B 2 ; s/ D ; Q B ; s D 1 K .B 2 ; s/ D .B 2 ; s/

2

(8.91)

where Z1 2

Q.B ; s/ D

Z1 d vv†f .v/

0

d v0 GQ vjv0 ; B 2 ; s '.B; v0 ; 0;/;

0

Z1 K.B ; s/ D

Z1 d vv†f .v/.1 ˇ/

2

0

(8.92)

d v0 GQ vjv0 ; B 2 ; s .v0 /;

(8.93)

0

and D.B 2 ; s/ D 1 K.B 2 ; s/ – sometimes referred to as the dispersion relation since its zeros are poles of .B Q 2 ; s/ that relate the discrete time decay constants ˛n 2 to B – has been introduced. It follows from this simple development that, if the 2 Q slowing-down equation, Eq. 8.88, can be solved for G.vjv 0 ; B ; s/, an expression 2 for the time-dependent response of a fission detector .B ; t/ in the experiment can be obtained via the inverse Laplace transform of Eq. 8.91. This was done using model cross sections and, first a synthetic (hydrogen-like) slowing-down kernel [55], and then the exact elastic slowing kernel [56, 57]. The key results, which were

8


415

essentially the same for the two studies, showed that the expressions for .B Q 2 ; s/ 2 2 had branch points at sbp D Œv†.v/ C vD.v/B min and poles sn .B / – at the zeros of D.B 2 ; s/ – and that as B 2 is increased – i.e., the assembly becomes smaller and farther subcritical – the last pole s0 .B 2 / “disappears” into the branch point. Prior to this there is a time-asymptotic exponentially decaying solution associated with the real, negative, isolated pole s0 .B 2 / D ˛0 .B 2 /. When that pole “disappears,” as B 2 is increased through a critical value B2 , it disappears only from the principal sheet of the Riemann surface for .B Q 2 ; s/. Actually, the pole bifurcates, as 2 2 B is increased through B , resulting in a complex conjugate pair of poles so˙ .B 2 / which were obtained by analytically continuing the expression for .B Q 2 ; s/ in the counter-clockwise (+) and clockwise (–) directions around the branch point sbp (see Fig. 8.8). The inverse Laplace transform is then given by a contour integral (path shown as dashed line in Fig. 8.8) that results from deforming the original integral along the Bromwich contour. When the integrand of this contour integral was expanded it led to time-dependence, for the detector response, of the form ˇ Cˇ ˇ t C 1 sC 2 t 2 ; Q B 2 ; t exp ˇsR 2 I

(8.94)

C on an intermediate time scale. Here, sR is the (negative) real part of the complex conjugate poles on the +1 and –1 Riemann sheets and sIC is the (very small) imaginary part.

s

Fig. 8.8 The Riemann surface with the trajectories of the bifurcated poles onto the adjacent (C1 and –1) sheets. The dashed curve represents the deformed contour of the Laplace inversion integral (Adapted from [54])

416

J. Dorning

If ˛.t/ is defined as ˛.t/ D .t/= .t/, P then on the intermediate time scale C C when Eq. 8.94 applies, it will be given by ˛.t/ jsR j .sIC /2 t, and sR will C 2 be the y-intercept of this local curve – an inclined plateau – and .sI / will be its slope. (See sketches of (a) log .B 2 ; t/ vs t and (b) ˛.t/ .B P 2 ; t /= .B 2 ; t / vs t in Fig. 8.9 (adapted from [54]).) Experimental data [58] comprising such an inclined plateau, indicating the existence of a pseudo mode and quasi-exponential die-away on a 10–40 s time scale in a 13 14 EURECA subcritical fast multiplying assembly, is shown in Fig. 8.10a (adapted from [54]). Analogous data [58] is shown in Fig. 8.10b (adapted from [54]), for a smaller, farther subcritical, 9 9 EURECA assembly, where the inclined plateau is less well-defined and much shorter-lived – 4–7 s [54]. As priorities in fast reactor development shifted toward other physics problems, e.g., sodium void effects in liquid metal fast breeder reactors (LMFBRs), and toward

a

b •

log rf

a (t) =

- rf rf

t

t

Fig. 8.9 (a) Sketch of the metastable (pseudomode) decay of a detector response. (b) Sketch of the time dependence of the logarithmic derivative, log .B 2 ; t/ (Adapted from [54])

a

0.20

b

a ( 1 ) msec

0.15 0.10 10

20

30

50 msec

40

0.9 0.8 0.7 0.6 0.5 1

2

3

4

5

6

7

8

9

10

11

12 msec

Fig. 8.10 (a) Experimental measurements of the “time-dependent” decay constant, ˛.t / .B P 2 ; t/= .B 2 ; t/ vs t for a EURECA 13 14 element assembly. Solid line is sketched in. (b) Similar data for a 9 9 EURECA assembly (Adapted from [54])

8


417

many pressing technological problems – and in general toward nuclear engineering problems directly related to LMFBR design in France, Japan, the Soviet Union and the United States, support waned for experiments of this type in bare fast subcritical assemblies; hence, the understanding of these and various other fundamental problems related to fast reactor physics problems did not progress nearly as far as the understanding that resulted from the earlier, more extensive studies of problems related to thermal reactor physics. Perhaps the budding resurgence of interest in fast reactors in France will change this, and put us on the path to developing a deep understanding of pulsed neutron experiments in subcritical fast reactor assemblies, and, therefore, a deep understanding of some important aspects of the kinetics of fast reactors.

8.5 Space–Time Reactor Kinetics In the early days of reactor development – the mid-1940s to the mid-1950s – when high speed, large digital computers did not operate at very high speeds, and did not have very large memory capacities, and were not very widely available – reactor physics analysis and reactor kinetics analysis was largely done using basic theoretical techniques and hand calculations. (In fact, when I was a graduate student, I met a computer programmer who told me his first job title had been “Computer.” Yes, he was a Computer! Early in his career he had done numerical calculations – by hand, and using a desk-top mechanical calculator – based on the expressions the physicists gave to him (e.g., the solution to the two-group, one-dimensional, steadystate diffusion equations in slab geometry).) Then, when the early antecedents of modern digital computers became available, e.g., in the US Navy’s nuclear reactor program, numerical solution techniques based on finite-difference schemes were programmed for the one- and two-group steady-state diffusion equations in slab geometry and for the one-group steady-state P-1 and P-3 equations in slab geometry. For some delightful reminiscences of that era the reader is strongly encouraged (dare I say, “commanded!”) to see the text of a wonderful after-dinner talk given by the late Dr. Ely Gelbard at an early American Nuclear Society Mathematics and Computational Division National Topical Meeting [59]. (I was there to enjoy it live, and it was one of the most interesting and entertaining talks I have ever heard!) As time passed, one-dimensional and later two-dimensional reactor kinetics codes based on finite-difference schemes in space and time were developed using the few-group and multigroup diffusion equations and the coupled delayed neutron precursor concentration equations. Some of these developments will be summarized briefly in Section 8.5.1. Because the memory capacity and speed of digital computers were still very limited in the late 1950s and early 1960s, alternative methods to those based on simple direct finite-difference schemes, were developed. These included methods based on variational functionals which were fairly widely used in reactor analysis (again, especially in the US Navy’s nuclear reactor program) in both steady-state

418

J. Dorning

reactor calculations and reactor kinetics calculations. Developed primarily during the early 1960s, they gradually evolved by the late 1960s into so-called synthesis methods in which the trial functions used in the variational formulations were the numerical outputs of the computer-generated solutions to related simpler problems. These methods were extensively developed and widely used in reactor analysis and design at both of the US Navy nuclear laboratories: Bettis Atomic Power Laboratory (BAPL) and Knolls Atomic Power Laboratory (KAPL). The basis for these methods will be reviewed in Section 8.5.2, and anomalies that arise in their application will be discussed briefly there. The early 1970s saw the development of another general class of computational methods for both reactor criticality calculations and reactor kinetics calculations. Also motivated by the goals of more accurate numerical solutions and more realistic representations of reactors using limited computer resources, they, to some extent, represented a logical compromise between the decomposition of a reactor into numerical cells (or boxes, or computational elements) characteristic of finite-difference methods and use of analytical solutions characteristic of variational methods. The basic idea of these so-called coarse mesh methods – some of which later bore the appellation nodal methods – was to decompose the reactor into computational elements or cells, as is done in finite-difference methods, and then rather than using the linearly truncated Taylor series approximation to the solution within the cell – as also is done in finite-difference schemes – introduce some more complicated, but better analytical approximations to the solution – as is done in variational methods. The distinction, of course, was that the solution is not approximated by global analytical functions over the whole reactor domain but rather by local analytical functions within each cell or computational element. This typically results in more accurate solutions – in comparison with those obtained using traditional finite difference methods – for a given cell size; or, conversely, it permits the use of larger, thus fewer, cells to achieve a solution of given accuracy requirements. Thus, specified solution accuracy was achieved using both less computer memory and less computer central processor unit (CPU) time on the mainframe computers of that era. Some of these coarse mesh [60] and nodal [61–67] methods for the multigroup diffusion equations, originally developed for criticality calculations, were soon afterwards extended to space–time reactor kinetics applications [62, 64, 65, 67] and some were also extended to the multigroup transport equations [68–73]. An early review of coarse-mesh and nodal methods for the diffusion equation appeared in 1979 [69], and an analogous review of nodal methods for the transport equation was written in 1985 [73]. Excellent, much more comprehensive, reviews of nodal methods for the diffusion equation and a related fuel-assembly homogenization procedure [74], and of nodal methods for the diffusion and transport equations [75] were published in 1986. The main ideas involved in the development of these coarse-mesh and nodal methods will be presented in Section 8.5.3 in which some results obtained using them will also be discussed. The success of these coarse mesh and nodal methods in achieving very high computational accuracy using very coarse meshes, or large computational elements or cells – in some cases as large as a light water reactor (LWR) fuel assembly

8


419

or larger – led to the need for a homogenization procedure that could be used to homogenize the fuel pins and associated lattice cells over a fuel assembly or a significant fraction thereof. One such homogenization procedure, called “equivalence theory,” in which an exact transport theory solution is introduced as a hypothetical solution on a computational element (fuel assembly) boundary, was introduced by K. Koebke [76], and subsequently used extensively and extended by G. Greenman, K. Smith and A. F. Henry [65] and K. Smith [74]. Others, based on asymptotic expansions were developed by several authors over a period of years. Most of these other homogenization theories began with an asymptotic expansion of the solution to the transport equation for a lattice of fuel cells in powers of the small parameter " D =R, where is an appropriate mean free path and R is a characteristic overall reactor dimension [77–81]. Clearly, this ratio is very small (of the order of 102 ) for a light water thermal reactor. One development, however, started from the asymptotic expansion of the solution to the diffusion equation in the small parameter © D L=R where L is the diffusion length; hence, this parameter is equally small (also of the order of 102 ) for a LWR. The essential ingredients that enter into these homogenization theories will be discussed in Section 8.5.4 in which a few results obtained using them in conjunction with nodal methods will be summarized. It is clear that, for different applications, it is appropriate to use different descriptions of reactor kinetics – the point reactor kinetics model, quasi-static models, space–time multigroup diffusion theory, and in some cases even space–time multigroup transport theory. It is also clear that, in some complicated transients, it is necessary to use a more precise model such as the time-dependent multigroup transport equations to adequately represent the portion of the transient in which, crucial local space–time variations of the neutron flux occur. However, it typically is very costly if such a precise computational representation is used throughout a long simulation of a transient. Thus, just as adaptive spatial grids and adaptive variable time steps are very effectively used in a whole host of computational physics and computational engineering problems – spatial grids in the neighborhood of shock-wave fronts in gas dynamics, refined time steps in phase transition problems, etc. – adaptive procedures to represent the evolution of the neutron flux during different epochs of a complex reactor transient by different models is an efficient way to proceed. Hence, this section on space–time reactor kinetics will close with a brief summary in Section 8.5.5 of two articles in which the development and application of precisely such an adaptive model procedure was reported [82, 83].

8.5.1 Finite-Difference Schemes for the Time-Dependent Multigroup Neutron Diffusion Equations The origins of numerical methods based on finite-difference schemes for the approximation of the solutions to differential equations go far back in time, as is evidenced by the names some of these schemes bear – for example, the forward Euler scheme and the backward Euler scheme, both of which are named after the

420

J. Dorning

Swiss mathematician Leonhard Euler, who lived from 1707 to 1783. These schemes were used extensively in “hand calculations” long before modern computers – or even desk-top mechanical calculators – had arrived on the scene. They are very simple and straightforward, mathematically well-understood, easy to program, and very widely applicable; and they are particularly suitable for the solution of diffusion equations. It is, therefore, not surprising that they were the first schemes employed in serious attempts at the numerical solution of reactor physics equations, particularly the neutron diffusion equation. In the 1950s and 1960s, and even now some 50 years later, most whole-core thermal reactor calculations begin, not from the multigroup transport equations, but rather from the few-group or multigroup diffusion equations 1 @

g .r; t / D r Dg .r; t /r g .r; t/ †a;g .r; t / g .r; t / vg @t C

G X

†s;g

g 0 .r; t/ g 0 .r; t /

g 0 D1 g 0 ¤g

C.1 ˇ/0;g

G X

†f;g 0 .r; t / g 0 .r; t/

g 0 D1

C

I X

i i;g Ci .r; t / C Qg .r; t /; g D 1; : : : ; G

(8.95)

i D1

and the coupled delayed neutron precursor equations G X @ Ci .r; t/ D ˇi i;g †f;g 0 g 0 .r; t / i Ci .r; t/ @t 0

i D 1; ; I ;

(8.96)

g D1

for reactor kinetics calculations, and from the steady-state version of these equations for reactor criticality calculations (k-calculations). Due to the dramatic difference between the time scale on which the prompt neutrons arrive – of the order of the prompt neutron lifetime – and that on which the delayed neutrons arrive – of the order of the precursor delay times – these equations are stiff differential equations. Thus, the time-eigenvalues associated with the delayed neutrons are much greater in magnitude than those associated with the prompt neutrons, and because of this an explicit finite difference scheme in time (such as the forward Euler scheme) would lead to an algorithm that would be numerically unstable in many applications. Hence, implicit schemes in time – such as the backward (implicit) Euler scheme (first order) or the Crank–Nicholson scheme (second order) – usually are used. Because most thermal reactor core layouts are very amenable to Cartesian geometry representation, the discussion that follows will use two-dimensional x–y

8


421

spatial coordinates. (The inclusion of the z-coordinate would be completely straightforward, but cumbrous.) Further, since there is no differential operator in the spatial variables, in the precursor equations, and their treatment is therefore straightforward, they will be omitted from the discussion. Finally, since the linear algebraic equations that result from applying difference schemes to the multigroup (and even the few group) diffusion equations comprise very large systems, which typically are solved via source iteration of the fission source, or the fission plus inscatter sources, the basic finite-difference scheme for the space–time solution is developed for just the g-th group equation. Hence, the essential part of Eqs. 8.95 and 8.96 that is relevant to this discussion is the time-dependent g-th group diffusion removal equation with a fixed source in x–y spatial coordinates @2 @2 1 @

g .x; y; t / D Dg .t/ 2 g .x; y; t / C Dg .t/ 2 g .x; y; t / vg @t @x @y (8.97) †r;g .t/ g .x; y; t / C Sg .x; y; t /; where the group removal cross section and diffusion coefficient have been taken to be uniform in space to avoid tedious and possibly confusing details; and Sg .x; y; t / represents the fission, inscatter, and precursor decay sources, along with any possible fixed sources. A simple central difference scheme applied to the second-order partial derivatives in x and y leads to 1 d D.t/ D.t/ D.t/

j:k .t/ D

j C1;k .t/ 2

j;k .t/ C

j 1;k .t/ 2 2 v dt x x x 2 D.t/ D.t/ D.t/ C

j;kC1 .t/ 2

j;k .t/ C

j;k1 .t/ 2 2 y y y 2 †r .t/ j;k .t/ C Sj;k .t/; j D 0; : : : ; J; k D 0; : : : ; K

(8.98)

where j;k .t/ D .xj ; yk ; t /, a uniform spatial mesh of cells xy has been used, and the group index g has been dropped. Clearing the 1=v and writing these equations in matrix form for the column vector .t/ of unknowns j;k .t/; j D 0; : : : ; J; k D 0; : : : ; K yields d

.t/ D L.t/ .t/ C S .t/; D dt

(8.99)

where the definitions of the square matrix L.t/ and the vector S .t/ follow from D Eq. 8.98 and are obvious. Finally, introducing a simple finite difference scheme for the time derivative leads to

` `1 D .1 /t` L ` ` t`1 L `1 `1 D

D

C.1 /t` S ` t`1 S `1 ; ` D 1; 2; 3; : : :

(8.100)

422

J. Dorning

where ` .t` /; L ` L.t` /; S ` S .t` /; t` is the discretized time, t` D

D

t` t`1 is the `th time step, and has been introduced so that a -weighted average of the right-hand-side of Eq. 8.99 appears in Eq. 8.100. Solving these equations for

` , the vector of fluxes at spatial grid points .xj ; yk / at time level t` in terms of

`1 and the sources S ` and S `1 yields 1

` D I .1 /t` L ` D D I t`1 L `1 `1 C .1 /t` S ` t`1 S `1 ; ` D 1; 2; 3; : : : D

D

(8.101) Here, the vector of sources S ` depends upon the group fluxes at the current (the `th) time level. Some of these group fluxes have been computed in the current outer source iteration and therefore are available; the others have to be approximated by their values at the previous time level t`1 . When , in this so-called theta-difference method in time, is set equal to unity, Eq. 8.101 becomes the result for the forward Euler scheme, which has a first-order global error in t (for a uniform time step) and is explicit (no matrix inversion required) – but often unstable. For equal to 0, this equation gives the timeadvancement operator for the backward Euler scheme, which also is first order, but implicit (matrix inversion required) and stable. Finally, for equal to one half, it gives the advancement operator for the Crank–Nicolson scheme, which has a second-order global error in t (for a uniform time step) and also is implicit and stable. Of course, there are many very important details – cell-centered versus celledge spatial finite-difference schemes, acceleration schemes such as coarse-mesh rebalancing, asymptotic source extrapolation, Wielandt (eigenvalue shift) iteration, etc. – that have been omitted here in this brief summary. However, many of these are discussed at some length in well-known textbooks [84–87]. The -difference scheme just summarized, albeit rather briefly, was the basis for the one-dimensional space–time reactor kinetics code WIGLE [88], developed at BAPL in the early 1960s – and also for TWIGLE [89], its extension later in the 1960s to two-dimensional space–time reactor kinetics. These codes were very widely and effectively used in the 1960s and 1970s; and their “guts” can still be found in some reactor kinetics and reactor dynamics software used today. Unfortunately, in order to achieve reasonable accuracy in calculations for LWRs – in which the neutron diffusion length is of the order of a couple of centimeters – very fine spatial meshes must be used when these standard finite difference schemes are employed. Thus, the vectors of group fluxes, ` in Eq. 8.101, become quite large even in two-dimensional calculations. For example, if only one hundred points were used in each direction (x and y) in a quarter-core calculation, that vector would have dimension 10,000 and the matrix that would have to be

8


423

inverted (in Eq. 8.101) in each group, in each iteration, at each time step would be 10; 000 10; 000. And for three-dimensional quarter-core calculations – which are important in order to properly represent the motion of control rods or blades in reactor kinetics calculations – these numbers become 2 106 and .2 106 / .2 106 /! Thus, in those days, when “large jobs” were run on (32K!) mainframe computers (e.g., the IBM 360) only at night or over weekends, two-dimensional LWR kinetics calculations were not done routinely! And three-dimensional calculations were not done at all – at least not using codes based on finite difference schemes. Rather, other methods that required less computer storage and time were developed both for criticality calculations and space–time reactor kinetics calculations. These included methods developed in the 1960s based on variational principles – in which the space-dependence of the neutron flux typically was expanded in functions defined over the whole reactor spatial domain – and later, in the 1970s and 1980s, methods that were capable of achieving very high numerical accuracy while using very large computational elements, or, equivalently, coarse meshes. The main ideas used in the development of these methods based on variational principles, along with some results obtained using them, will be summarized in the next section. And the key steps in the formulation of the so-called coarse-mesh and nodal methods and some examples of their application will be described in the section after that.

8.5.2 Variational, Modal, Synthesis, and Related Methods for the Time-Dependent Multigroup Diffusion Equations The development of variational methods, (variational) synthesis methods and socalled modal methods for space–time reactor kinetics (and reactor statics) begins from a simple variational principle [90–94]; hence, it seems apropos to include here a brief review of the use of these principles to develop techniques for obtaining approximate solutions to differential equations, integral equations, and other functional equations. Even though none of the equations of interest in space–time reactor kinetics – except the one-group steady-state diffusion equation – is selfadjoint it is convenient for expository purposes to begin from a linear self-adjoint equation L D S

(8.102)

and introduce the classic functional [90–92] GŒ' D .'; L'/ 2.'; S/:

(8.103)

Here, all the variables are taken to be real, and the real inner product is given by the integral over the multidimensional domain X over which the unknown function '.x/ is defined. The function that minimizes this functional is the '.x/ that is the

424

J. Dorning

solution to Eq. 8.102. This becomes obvious when the variation of GŒ' with respect to ' is set equal to zero ıGŒ' D .ı'; L'/ C .'; Lı'/ 2.ı'; S/ D 0;

(8.104)

D .ı'; L'/ C .L'; ı'/ 2.ı'; S/ D 0; D 2.ı'; fL' S g/ D 0;

(8.105) (8.106)

Here, ı' is the variation in ', and since it is arbitrary, the inner product is zero only if '.x/ is the solution to Eq. 8.102. (If this development began from Eq. 8.103 in which GŒ' represented some physical quantity – for example, the energy in a mechanical system – its minimization would yield Eq. 8.102 as the equation(s) of motion for the system, the so-called Euler-LaGrange equation(s) corresponding to the functional GŒ'. The second variation of GŒ' with respect to ' is 2.ı'; Lı'/ where L' S D 0 has been used; hence, if L is a positive operator the extremum that follows from ıGŒ' D 0 is a minimum. In practice, when the original operator in Eq. 8.102 is known – as is the case in reactor kinetics and reactor physics in general – an approximation to its solution can be generated by minimizing GŒ', not with respect to an arbitrary variation of ', but rather with respect to the variation of a trial function – comprising, for example, a linear combination of acceptable functions, i.e., functions that belong to some specific class. This procedure is made very explicit by introducing the trial function 'T .x/ D

N X

Cn n .x/;

(8.107)

nD1

into the functional GŒ' and then minimizing this functional by setting all its partial derivatives with respect to the Cn equal to zero – or equivalently setting its first variation fıGŒ T =ıCn gıCn ; n D 1; : : : ; N with respect to all the Cn equal to zero – which leads to N X

.n ; Lm /Cm .n ; S / D 0; n D 1; : : : ; N;

(8.108)

mD1

where a common factor of 2 has been cancelled, and the facts that L is self-adjoint and the inner product is real have been used. This is a simple nonhomogeneous linear system of algebraic equations which can be inverted as long as it does not have a zero eigenvalue, to obtain the Cn ; n D 1; : : : ; N and therefore, via Eq. 8.107, the specific 'T .x/ that minimizes the functional GŒ'T over functions of the form given by that equation, and thereby provides an approximation of this form to the solution of Eq. 8.102. This procedure is known as the Rayleigh–Ritz method [91], and it is equivalent to the Bubnov–Galerkin method [91] of approximation to solutions in linear function spaces. In that method, the solution to Eq. 8.102 is simply approximated by the expansion '.x/ Š 'N .x/

N X nD1

Cn n .x/;

(8.109)

8


425

in the set of functions n .x/; n D 1; : : : ; N , which is substituted directly into Eq. 8.102; then, the inner product of the resulting equation with each of the n .x/; n D 1; : : : ; N , is separately set equal to zero, forcing the residual RN .x/ D L'N .x/ S.x/ to be orthogonal to each of the expansion functions n .x/; n D 1; : : : ; N . The final equations, therefore, are identical to Eq. 8.108; hence, this method is sometimes called the Ritz–Galerkin method – even though the lives of the two mathematicians did not overlap. (Equation 8.108 is equivalent to setting the gradient of GŒ'T with respect to the vector C of coefficients Cn ; n D 1; : : : ; N , equal to zero, i.e., rC GŒ'n D 0.) This, of course, establishes that GŒ'T undergoes an extremum for the resulting value of C . To insure that this extremum actually is a minimum – and not a maximum or inflection point – the Hessian rC rC GŒ'T also should be calculated and evaluated at the value of C that corresponds to the extremum to show that it is positive and, therefore that GŒ'T has been minimized, not maximized. But this long and tedious step is usually omitted in real-world calculations. Now that the development for self-adjoint operators has been summarized, the extension to non-self-adjoint operators can be made clear very easily. This development – which applies to the time-dependent (and steady-state) multigroup diffusion equations and multigroup transport equations, and the time- and energy-dependent diffusion and transport equations as well – begins from the linear non-self-adjoint operator equation A' D S; (8.110) and the related adjoint equation A ' D S ;

(8.111)

where, for example, A might be a nonsymmetric matrix of non-self-adjoint integrodifferential operators, and A , of course, is its adjoint. As for the self-adjoint case a functional is introduced and its variation is set equal to zero. Here, however, the classic functional GŒ', given by Eq. 8.103, is replaced by its generalization to the non-self-adjoint case – the Roussopoulos functional [93, 94] (Reference [94], cited here, is a somewhat obscure, but lovely, introduction to variational methods in general and variational methods for reactor calculations in particular; and every one interested in reactor calculations and calculational methods should have a copy in her or his library. (My copy is a not-very-good photocopy of an old mimeographed copy; nevertheless, I cherish it!)) i h F ' ; ' D ' ; A' ' ; S S ; ' ;

(8.112)

where, in general, all the variables and the operator A may be complex and the inner product is the standard complex inner product over the multidimensional domain X. The variation of F Œ' ; ' with respect to both ' and ' , set equal to zero, is h i ıF ' ; ' D ı' ; fA' S g C fA ' S g; ı' D 0;

(8.113)

426

J. Dorning

and since both ı' and ı' are arbitrary variations this leads to Eqs. 8.110 and 8.111 for the so-called forward and adjoint solutions, respectively, as the generalized Euler–LaGrange equations. (Analogously to the self-adjoint case, if the original equations for ' and ' , Eqs. 8.110 and 8.111, had not been known and this development began from the functional F Œ' ; ', setting its first variation equal to zero would have generated these equations.) Now, approximations to the solutions to the forward equation and the adjoint equation can be generated via steps closely analogous to those just summarized for the self-adjoint case. Introducing trial functions for both the forward and adjoint solutions 'T .x/ D

N X

Cn n .x/;

(8.114)

Ck k .x/;

(8.115)

nD1

and 'T .x/ D

N X kD1

into the Roussopoulos functional and setting its first variation with respect to all the Ck and all the Cn equal to zero N n X

N n h o h o i i X @F ' ; ' =@Ck ıCk C @F ' ; ' =@Cn ıCn D 0;

(8.116)

nD1

kD1

or equivalently, since all the ıCk and all the ıCn are arbitrary, setting all the par tial derivatives of F Œ' ; ' with respect to both the Ck ; k D 1; N and the Cn ; n D 1; N equal to zero yields N X

` ; An C

n

` ; S D 0; ` D 1; : : : ; N ;

(8.117)

nD1

and

N X k ; Am C m S ; m D 0; m D 1; : : : ; N :

(8.118)

kD1

Equations 8.117 form a linear nonhomogeneous algebraic system for the C` ; ` D 1; : : : ; N , the solution of which gives the forward trial function 'T .x/ in Eq. 8.114 which, along with the adjoint trial function generated via Eq. 8.115 for the Cm ; m D 1; : : : ; N from the solution to Eq. 8.118, renders the functional F Œ' ; ' stationary and, therefore provides an approximation to the solution of the original equation, Eq. 8.110. And, naturally, if an approximation 'T .x/ to the solution of

the adjoint equation is desired, Eq. 8.118 can be solved for C m ; m D 1; : : : ; N . Here, in the non-self-adjoint case, since the functional is rendered stationary with

8


427

respect to both ı' and ı' it is not known whether the extremum in the functional that corresponds to the approximate solutions for '.x/ and ' .x/ is a minimum, a maximum or a saddle or inflection point. Just as the Rayleigh–Ritz equations, Eqs. 8.108, that followed from the classic variational functional for self-adjoint operators, can also be derived via the simple Bubnov–Galerkin method, so also can Eqs. 8.117 and 8.118 be developed via simple linear function space approximation methods. If the solution to Eqs. 8.110 and 8.111 are approximated by the expansions '.x/ Š 'N .x/

N X

Cn n .x/;

(8.119)

Ck k .x/;

(8.120)

nD1

and ' .x/ Š N .x/

N X kD1

and this 'N .x/ is substituted into Eq. 8.110, and the resulting residual RN .x/ A'N .x/ S.x/ is required to be orthogonal to each of the adjoint expansion func tions 'k .x/; k D 1; : : : ; N , Eqs. 8.117 immediately follow. Similarly, if 'N .x/ is .x/ substituted into Eq. 8.111, and the residual thus generated R .x/ D A 'N S .x/ is made orthogonal to each of the forward expansion functions n .x/; n D 1; : : : ; N , Eqs. 8.118 result immediately. Further, if the residual R.x/ is required to be orthogonal to each of a more general set of so-called weight functions !k .x/; k D 1; : : : ; N , the result is Eqs. 8.117 with the ` .x/ replaced by the !` .x/. These equations correspond to the weighted residuals method to generate an approximation to the solution of Eq. 8.110. Clearly, if the !k .x/ are taken as the expansion functions k .x/ themselves, the result is the Bubnov–Galerkin method (mentioned above) applied to the non-self-adjoint equation, Eq. 8.110. Finally, if the weight functions !k .x/ are taken to be the residuals themselves, the result is known as the least squares method. These basic ideas on the development of variational methods for non-self-adjoint problems can be applied directly to space–time reactor kinetics calculations and, of course, also to reactor criticality calculations. To facilitate this, the operator A in Eq. 8.110 and the inner product in Eq. 8.112 were kept quite general. For example, A might represent the operator in the time-dependent multigroup diffusion equations and coupled precursor equations, Eqs. 8.95 and 8.96, including the time derivatives – in which case the inner product would include not only the summation over all g D 1; ; G and all i D 1; : : : ; I and an integration over all r in the reactor volume V , but also an integral over time from the initial time t0 to the final time tf . Or, A could represent only the operator on the right-hand sides of these equations, thereby excluding the time derivatives – in which case the integral over time would not be included in the inner product. When variational methods were first introduced into reactor statics analysis and reactor kinetics analysis, analytical expressions were used for the terms in the sums

428

J. Dorning

that appeared as the trial functions. However, it was soon realized that more accurate solutions to practical reactor physics and reactor kinetics problems would result if the trial functions were comprised of functions that were solutions to very closely related, realistic, model problems. These solutions, of course, could be provided by generating numerical solutions to problems that are based on an essentially identical representation of the reactor but which are cast in a simpler context. For example, if the reactor kinetics problem that was to be solved were a control-rod-driven transient in a thermal reactor modeled by a two-group, three-dimensional, space-time, diffusion-theory description in a specific geometric representation, the trial function terms to be used in the context of the transient calculation based on the variational principle might be those generated in three reactor statics (criticality) calculations using the same two-group, three-dimensional representation of the reactor: one with the rods out, one with the rods half-way in and one with the rods all the way in. When the trial functions used in solutions based on variational principles were comprised of numerical solutions, such as in this example, the methods were called “variational synthesis methods,” or more commonly, simply synthesis methods [15]. Since numerical trial functions used in these synthesis methods do not, in general, have the smoothness properties (continuity, continuous first derivatives, etc.) required in the formal development of the underlying variational (or weighted residuals) methods, these trial functions were not admissible trial functions, and the solutions obtained using these synthesis methods sometimes exhibited anomalous behavior, such as unphysical spikes or other discontinuities, for example, in the calculated time evolution of the neutron flux in a reactor kinetics computation. However, mathematically literate reactor engineers were well aware of all this, and they usually simply ignored these unphysical anomalies. In the lexicon of variational methods for reactor analysis, the class of methods known as modal methods is simply a special case in which the trial functions are comprised of “modes” or eigenfuctions that arise in related problems. In fact, the example given above as a synthesis method could also be described as a modal method, since the trial function there was comprised of the fundamental (numerical) modes or eigenfunctions from three criticality calculations. In practice, it often is very useful to use different functions in the trial functions in different spatial regions of the reactor (or at different times in a transient). In such cases the functionals given in Eqs. 8.103 and 8.112 are not adequate because they cannot accommodate discontinuous trial functions. However, these functionals – most importantly, since it is relevant to non-self-adjoint operators, the Roussopoulos functional given in Eq. 8.112 – can easily be extended to include discontinuous trial functions (as well as trial functions that do not satisfy the (external) boundary conditions or the initial (or final) conditions). Although the development of these more general functionals is fairly straightforward, the final expressions are very long and very complicated [95]; hence, they are not included here. These methods generally bore the name “multichannel synthesis methods,” since different (numerical) trial functions were used in different regions or “channels” of the reactor and they are still used in statics and kinetics calculations done in support of thermal reactor design today – in the twenty-first century. (For a comparison of the results of

8


429

a sequence of multichannel synthesis calculations, done using successively larger numbers of channels, with the results of a fine-mesh finite-difference calculation, for a transient in a simple reflected slab reactor model, see [15], Figs. 3.1 and 3.2 and the related discussion on pages 77–80.)

8.5.3 Coarse-Mesh and Nodal Methods for Space–Time Reactor Kinetics The quest to develop advanced computational methods that resulted in better accuracy – in the numerical solutions generated in support of reactor statics and kinetics analysis – and/or that utilized less computer resources for given accuracy requirements, first led from traditional finite difference methods to methods based on variational principles such as multichannel synthesis methods. Each of these two classes of methods had its advantages; but each also had its disadvantages. The finite-difference methods required very large numbers of discrete-variable unknowns for accurate multidimensional calculations, resulting in large computer storage requirements and long (cpu) running times; and the variational synthesis methods had the problem of the occasional unphysical discontinuities when practical numerical trial functions were used, and, as the number of channels used in the multichannel synthesis methods was increased, the number of discrete-variable unknowns increased – resulting in increasing demands on computer resources. Thus, the need for other, more efficient, computational methods continued to exist. In response to this continued need, a new class of numerical methods described, in general, as “coarse-mesh methods” was developed. The common feature of these methods was the use of large computational elements – large in comparison with the cells used in traditional first-order finite-difference methods – combined with various techniques that would still lead to numerical results on the associated “coarse meshes” that were as accurate as, or more accurate than those achieved by finitedifference methods using small cells, or “fine meshes.” For example, in one such method an analytical solution (a local Green’s function) was used within a coarsemesh two- or three-dimensional computational element to give the outgoing partial current – on a surface of the element – that results from the incoming partial currents on all the surfaces of the element and the sources within the element [60]. Other coarse-mesh methods used higher-order numerical approximations within a coarse-mesh computational element to achieve the same objective. However, in a sense, some of these methods were “too accurate” in that the numerical solutions to which they led had much, much smaller errors than the errors in the nuclear cross sections and other nuclear data in the multigroup diffusion equations that were being solved. Thus, there was a need for coarse-mesh methods that could be less accurate (but, of course, still very accurate) but would require the computation of fewer discrete-variable unknowns per computational element. This need was answered by the appearance of a class of coarse-mesh methods to which we shall refer to as “transverse-integration nodal methods.”

430

J. Dorning

Since the high computational efficiencies of these methods are derived from the fact that they result in high accuracy while using large computational elements in the two- or three-dimensional spatial domain, their benefits are essentially the same in the contexts of reactor kinetics calculations and reactor statics (criticality) calculations. Moreover, because in their applications to reactor kinetics they are usually combined with implicit finite-difference methods in time, the formal development of these methods is essentially the same for the two applications. (This is not the case in the implementation of these ideas to computational heat transfer [96] and computational fluid dynamics [97, 98] – where the partial differential equations are not normally stiff in the time variable, as they are in reactor kinetics due to the dramatically different time scales associated with the prompt neutron lifetime and the delayed neutron precursor decay times.) Hence, to illustrate the basic ideas of the development of transverse-integrated nodal methods for the multigroup, multidimensional, space–time reactor kinetics equations, it is sufficient to illustrate them in the context of the steady-state diffusion equation. Further, since the multigroup neutron diffusion equations are solved group-by-group using power iterations [85] in nodal methods applications to reactor kinetics, only the one-group diffusion equation need be considered. Finally, all the key ideas can be made clear in the setting of a two spatial dimension model. Hence, we begin the development here from the two-dimensional one-group diffusion equation at time tk with spatially uniform cross sections and diffusion coefficient within a rectangular computational element (or “node”) Dg

@2 @2 ' .x; y; t / D 'g .x; y; tk / C †rg 'g .x; y; tk / g g k @x 2 @y 2 D Sg .x; y; tk /; x 2 .ai ; Cai /; y 2 .bj ; Cbj /;

(8.121)

i D 1; : : : ; I; j D 1; : : : ; J; g D 1; : : : ; G; k D 1; 2; : : : where x and y are local spatial coordinates in the ij-th node which is of width .ai ; Cai / in the horizontal direction and height .bj ; Cbj / in the vertical direction, Dg is the group diffusion coefficient, 'g .x; y; tk / is the group scalar flux, Sg .x; y; tk / is the source of neutrons in the g-th group, g is the group index, and G is the number of groups. In the context of reactor kinetics, †rg is the g-th group removal cross section plus 1=.vg tk / where vg is the g-th group neutron velocity and tk D tk tk1 is the k-th time step in the backward difference scheme for the time derivative, and Sg .x; y; tk / includes the g-th group inscatter and fission sources (the latter of which explicitly gives the prompt neutron fission source and the precursordependent delayed neutron fission source) plus Œ1=.vg tk /'g .x; y; tk1 / which is the other term from the backward difference scheme for the time derivative and which depends upon the previously computed neutron flux at tk1 . Suppressing the group subscript g and the time argument tk the starting point becomes D

@2 @2 '.x; y/ D '.x; y/ C †r '.x; y/ D S.x; y/: @x 2 @y 2

(8.122)

8


431

To develop a transverse-integrated (or equivalently, a transverse-averaged) nodal method, this partial differential equation is first integrated over one of the independent spatial variables within the ij -th computational element, e.g., y 2 .bj ; Cbj / to obtain an ODE in the other variable, e.g., x 2 .ai ; Cai / D

d2 y 1 y ŒJy .x; Cb/ Jy .x; b/ ' .x/ C †r ' y .x/ D S .x/ 2b dx2 y S .x/; (8.123)

where the computational element (or “node”) index j has been suppressed on b, and the “barred” (transverse-averaged) quantities (neutron flux ' y .x/ and source y S .x/) are defined within the node by 1 Q .x/ D 2b y

Z

Cb

Q.x; y/dy; Q D '; S;

(8.124)

b

and Jy .x; ˙b/ are the net neutron currents at the ij -node surfaces y D ˙b. Then Eq. 8.122 is separately transverse-integrated (or transverse-averaged) over the other spatial independent variable within the ij -th node, x 2 .ai ; Cai / to obtain an ODE in the second independent variable y 2 .bj ; Cbj / D

1 d2 x x ŒJx .Ca; y/ Jx .a; y/ ' .y/ C †r ' x .y/ D S .y/ dy 2 2a x

S .y/;

(8.125)

where the node index i has been suppressed on a and the “barred” quantities are defined within the ij -th node via definitions analogous to those in Eq. 8.124 x

Q .y/ D

1 2a

Z

Ca

Q.x; y/dx; Q D '; S;

(8.126)

a

and, of course, Jx .˙a; y/ are the net currents at the ij -node surfaces x D ˙a. Each of these two transverse-averaged ODEs, Eqs. 8.123 and 8.125 can easily be solved to determine the transverse-averaged flux, or net current or partial current, at a node surface in terms of the same quantity on the opposite surface and the source within the node. y Since the generalized source term S .x/ in Eq. 8.123 depends upon the net currents on the transverse node surfaces, Jy .x; Cb/ and Jy .x; b/, the solution to y this equation, e.g., for the y-average of the net current at x D Ca; J .Ca/ is only a formal solution. However, these transverse surface currents, Jy .x; Cb/ and Jy .x; b/, can be approximated by the solutions to Eq. 8.125 in the two vertix x cally adjacent nodes (in the y-direction) for J y .Cb/ and J y .b/. Moreover, this “constant” approximation for the x-dependence of these transverse surface currents, which is quite crude, can be replaced by a “quadratic transverse surface fit” in which

432

J. Dorning x

the solutions for J y .Cb/ to Eq. 8.125 on the upper surface of three horizontally adjacent nodes, the I –1, j node, the i; j node and the i C 1, j node are used on the interval x 2 .3a; C3a/ to develop a quadratic approximation for Jy .x; Cb/ x that correctly preserves the averages J y .Cb/ over the upper surfaces of each of these three nodes. Thus, when a similar approximation is introduced for Jy .x; b/ in Eq. 8.123, and analogous transverse quadratic fits are used for Jx .Ca; y/ and Jx .a; y/ in Eq. 8.125, the analytical solutions to Eqs. 8.123 and 8.125 can be used y to generate a set of linear algebraic equations for the node surface quantities J x .˙a/ x and J y .˙b/ on all the nodes. These equations, however, are underdetermined, and to increase the number of equations to equal the number of unknowns, two physically obvious and mathematically important constraints must be imposed on each x;y must be the computational element or node: (1) the double average of the flux ' same – whether it is calculated from ' y .x/ or ' x .y/, and (2) the double average x;y

of the generalized source S also must be the same – independently of whether y x it is computed from S .x/ or S .y/ [99]. Here the double averaged quantities are defined by Z Ca Z Cb x;y 1 1 y x Q D Q .x/dx D Q .y/dy; Q D '; S: (8.127) 2a a 2b b When these so-called constraint equations are added the linear system for the node surface-averaged currents becomes well-posed and can be solved iteratively within each energy group g and then for all the energy groups G at each successive time step tk . Various transverse-averaged (or transverse-integrated) nodal methods – all of which are capable of very high accuracies on coarse (spatial) computational meshes – were developed, most in the 1970s and 1980s, based on various solutions to Eqs. 8.123 and 8.125. The “nodal expansion method” [66, 67] was based on solutions to these equations obtained using approximate expansions, and the “analytic method” [61, 65, 74] was developed, with thermal reactor applications in mind, by solving Eq. 8.123 analytically for two simultaneous energy groups – and, of course, also solving Eq. 8.125 for two simultaneous energy groups. In the “nodal Green’s function method,” [63, 64] the net currents in the generalized source terms on the right-hand sides of Eqs. 8.123 and 8.125 were rewritten in terms of the partial currents in opposite directions on the same surfaces, and Eqs. 8.123 and 8.125 were solved analytically for a single group (in anticipation of general multigroup applications). (Green’s functions were used for historical reasons to obtain the very simple analytical solutions in that development because they had been used earlier in the derivation of the partial current balance method in which their use was convenient in the construction of the analytical solutions to two- and three-dimensional diffusion equations in computational elements – since transverse averages were not introduced in the formulation of that method [60]. Because the nodal Green’s function method used the partial currents on the node surfaces as the discrete-variable unknowns, the “hybrid nodal Green’s function method” [100] – which employed the net

8


433

currents on these surfaces, as in Eqs. 8.123 and 8.125 – was subsequently developed, reducing the number of unknowns that had to be computed by a factor of two [100]. Later, a higher-order nodal method – in which transverseaveraged group diffusion equations for higher-order transverse moments of the flux are introduced [69] – was developed as a significant extension of this method [101]. All these transverse-averaged nodal methods derive their high computational efficiency from the high numerical accuracy they achieve on coarse spatial meshes in reactor statics calculations. Coarse-mesh methods in general, and nodal methods in particular, yield highly accurate results for LWR calculations using computational elements as large as a fuel assembly – 10–20 cm on each side in two-dimensional models. Thus, in comparison with traditional finite-difference methods, which require the use of a least one computational element or “cell” for each fuel-pin lattice cell – one cm on each side – these methods require the solution for about one hundred times fewer discrete-variable unknowns, and thereby also require about one hundred times less computing time, and a great deal less computer storage, to achieve comparable accuracy in a two-dimensional calculation. And in three-dimensional calculations in LWR applications, this factor increases from approximately 102 to 103 . Thus the decrease in the number of unknowns, and the corresponding decrease in the required computing time increases from a factor of about 100 to a factor of about 1,000. To put this reduction in computing time, for comparable practical accuracy requirements, in perspective, it is useful to recognize that a factor of 1,000 corresponds to the replacement of a calculation that takes a whole day (24 h) by one that takes less than 1 1=2 min! (That is, it is equivalent to the difference between a “work day and overnight sleep at home” and a “very, very short coffee break!”) Because a spatial solution calculation (two-group, few-group, or multigroup) must be carried out at the end of each time step in a reactor kinetics computation, the high computational efficiencies of coarse-mesh and nodal methods in reactor statics calculations carries over directly to reactor kinetics applications. Numerous examples that illustrate these advantages in such applications are available in the articles that reported the original development of these methods for reactor statics and in later articles in which applications of these methods to reactor kinetics computations are reported [60–67]. The success of these coarse-mesh and nodal methods in achieving high accuracy on very coarse meshes was dramatically underscored by the fact that the errors that resulted in the calculations done using these methods were typically very much smaller than those that resulted from the procedures used to generate the nuclear data (i.e., cross sections and diffusion coefficients) homogenized over the coarse computational elements. Thus, the largest component of the error in the final reactor statics calculation or reactor kinetics computation was due, not to the nodal methods calculation, but rather to the homogenization procedure. Not surprisingly, this circumstance led to efforts to develop more accurate homogenization methods in general and more accurate homogenization methods for use in conjunction with coarse-mesh and nodal computational methods in particular. In the next section,

434

J. Dorning

a brief summary of two such homogenization procedures will be given in order to complete this summary of coarse-mesh and nodal methods for computations in space–time reactor kinetics analysis. Before closing this section on coarse-mesh and nodal methods for the timedependent multigroup diffusion equations, however, it should be mentioned that transverse-integration nodal methods were extended from the neutron diffusion equation to the neutron transport equation in the 1970s and 1980s [68–73, 75], and higher-order nodal methods also were developed for the transport equation [102]. Discussions of nodal methods for solving the transport equation, as well as other methods for the transport equation based on the use of local (within cell) analytical solutions, are presented in Section 1.4 of this book.

8.5.4 Homogenization Theories for Space–Time Kinetics Calculations and for Point Kinetics Calculations The high computational accuracies on coarse computational meshes that resulted from transverse-averaged nodal methods, and other coarse-mesh methods, for the steady-state and time-dependent multigroup, multidimensional neutron diffusion equations made it very obvious that better homogenization procedures would have to be developed to provide more accurate cross sections and diffusion coefficients, homogenized over the large computational elements, if the full advantage of these advanced computational methods were to be realized. Two different “homogenization theories” were developed in response to this need. Although they are based on fundamentally different approaches, the final homogenization procedures have some features in common. These two “approaches to homogenization” – which seems a more appropriate description than “homogenization theory,” which sounds a little grandiose, but which, nevertheless, has found its way into the literature – will be summarized here because they have been used in conjunction with transverseaveraged nodal methods for the neutron diffusion equation. In the first of these approaches, called exact equivalence theory, the exact solution (diffusion theory, transport theory, etc.) for a fully heterogeneous coarse-mesh computational element – usually a fuel assembly or a quarter of a fuel assembly in a two-dimensional calculation – is introduced. This exact “heterogeneous” solution is then used to generate node-homogenized (coarse computational element homogenized) cross sections and “direction-dependent” diffusion coefficients, and assembly “heterogeneity factors” for use in the homogenized diffusion equation so that it will yield a solution for keff and the node-averaged scalar fluxes that are equivalent to those which would result from the exact equations for the heterogeneous reactor [74, 76, 103]. Hence the name “exact equivalence theory.” The introduction of direction-dependent diffusion coefficients and an assembly heterogeneity factor for each node provides the additional degrees of freedom that make it possible for the homogenized diffusion equation to replicate keff and the node-average scalar fluxes that would result from the exact equations for the full heterogeneous fuel assembly or quarter fuel assembly.

8


435

Of course in a practical application, the exact solution for the heterogeneous node in the reactor being modeled is not available. (If it were, there would be no need to do the diffusion theory calculation!) Hence, some other “reference” heterogeneous solution must be substituted for it. Typically, the replacement solution is one obtained for the heterogeneous computational element surrounded by an infinite lattice of replicas of that element – which in practice requires only a criticality calculation for the isolated single heterogeneous element with zero net current boundary conditions – an infinite medium k-calculation for the heterogeneous node. This, of course, is a very reasonable practical substitution; however, it clearly is not the exact solution for the heterogeneous computational element in the actual environment of the original heterogeneous reactor model developed for the statics or kinetics calculation. Of course, better approximations for the solution for the scalar flux within the heterogeneous node can be generated by surrounding the heterogeneous node by a buffer zone of nodes from the heterogeneous reactor; however, the result would still not be the solution for the heterogeneous node in the original heterogeneous reactor. This exact equivalence theory was modified and extended to develop an approach to homogenization that came to be known as generalized equivalence theory [74, 104]. A reference solution for the flux in the coarse computational element also is used in the homogenization procedure that results from this approach; and, in practice, it usually is the solution that is obtained from a calculation of the type just described. However, instead of a single flux heterogeneity factor for each node (for each direction) – equal on opposite sides of the node, two flux discontinuity factors – or flux jump conditions – on the two opposite surfaces of the node, are introduced to provide the degrees of freedom that make it possible for keff and the node-average scalar fluxes to be preserved by the node-homogenized diffusion equation. The resulting homogenized cross sections and diffusion tensor are given by their space-dependent values in the interior of the heterogeneous node weighted by the node-interior scalar flux and integrated over the node, and of course, divided by the integral of this flux over the node; and the flux discontinuity factors are used in the whole-reactor node-homogenized diffusion theory calculation to provide prescribed interface jump conditions on the surfaces between adjacent nodes. After that calculation is completed, an approximation of the node-interior space-dependent scalar flux can be synthesized by multiplying the normalized heterogeneous flux for each node by the node-averaged flux computed in the diffusion calculation. The second of the two general approaches to homogenization for coarse-mesh and nodal diffusion methods that will be summarized here is based on the asymptotic expansion of the solution for the angular flux of the transport equation combined with the use of multiple spatial scales. The basic idea of the expansion of the angular flux in the small parameter, defined by the ratio of the particle mean free path to the characteristic dimension of the macroscopic physical system, goes back quite far in history. In fact, it can be found in the Hilbert expansion solution to the nonlinear Boltzmann equation in the kinetic theory of molecular gases, where Hilbert used such an asymptotic expansion to show that leading-order solutions to the Boltzmann equation can be obtained in the small mean free path limit by solving hydrodynamics equations [105]. An expansion of this type, in the ratio of the

436

J. Dorning

neutron scattering mean free path to the characteristic reactor dimension, was used to develop a fuel-pin lattice cell homogenization procedure [77, 78]. Starting from the neutron transport equation for a reactor comprising a spatially periodic assemblage of identical lattice cells, such an expansion of the angular flux was introduced, along with two spatial scales – one corresponding to the dimension of the lattice cell (about one cm for a LWR, and therefore about one mean free path for such a reactor), and the other corresponding to the overall reactor dimension (which is on the order of meters). The solvability conditions, due to the Fredholm alternative applied to the hierarchy of infinite-medium spatially periodic lattice-cell transport equations that result, give rise to the leading-order solution as the product of the transport theory solution for the lattice-cell and the solution to the resulting lattice-cell-homogenized diffusion equation for the whole reactor. This development was extended in three ways for use with coarse-mesh and nodal diffusion methods. The first, and most straightforward, was simply based on the replacement of the fuel-pin lattice cells by heterogeneous fuel assemblies – and relaxing the perfect periodicity assumption to near-periodicity so that slow variations in space, such as those that result from nonuniform burn-up, could be accommodated [79, 80]. The result, of course, was a fuel-assembly homogenized diffusion equation to be used in the coarse-mesh or nodal diffusion calculation for the whole reactor. In addition the continuity – of the leading-order term in the asymptotic expansion of the transport theory solution – imposed there led to flux discontinuity conditions at the node interfaces in the fuel-assembly-homogenized diffusion equation. These, of course are reminiscent of the flux discontinuity factors proposed in generalized equivalence theory (see above). Conversely, the expressions for the assembly-homogenized cross sections and diffusion coefficient are not given by the transport-flux-weighted averages of these space-dependent quantities over an assembly; rather, they are weighted by both the transport theory flux and the transport theory adjoint flux (or importance) – a result that follows directly from the solvability conditions mentioned above. The second of the three ways the lattice-cell homogenization development was extended, in some ways seems to be a step backward. Instead of beginning from a transport theory description of the whole reactor, it started from a diffusion theory description of the whole heterogeneous system and the expansion was in powers of the ratio of the diffusion length – also of the order of one cm for a LWR – to the characteristic reactor dimension [106]. This was done because, in actuality, the initial model for whole LWR statics or kinetics calculations often is not the transport theory description of the heterogeneous system, but rather, the diffusion theory description. The results of this two-spatial-scale diffusion-theory-based development of procedures for fuel-assembly homogenization are completely analogous to those summarized in the preceding paragraph for the two-spatial-scale transport-theorybased development – except, naturally, the transport theory solutions and adjoint solutions for the heterogeneous fuel assemblies are replaced by the analogous diffusion theory solutions and adjoint solutions in the expressions for the homogenized cross sections and diffusion coefficients and for the flux discontinuity factors.

8


437

The third, and most important, extension of the lattice-cell homogenization development [81], used the asymptotic expansion of the transport theory solution for the full heterogeneous reactor in the ratio of the scattering mean free path to the characteristic reactor dimension. However, three, rather than two, spatial scales were introduced. The first was the lattice-cell scale – of the order of one cm for a LWR; the second was the fuel-assembly scale – of the order of 10 cm; and the third was the whole reactor scale – of the order of 100 cm. This made possible a single analytical development that led simultaneously to self-consistent procedures for the homogenization of heterogeneous fuel-pin lattice cells and the homogenization of heterogeneous fuel assemblies. The results were lattice-cell-homogenized diffusion equations for each fuel assembly, and a fuel-assembly-homogenized diffusion equation for the whole reactor. The lattice-cell homogenized fuel-assembly diffusion equations have cross sections and diffusion coefficients that are given for each cell by the cell averages of the cross sections and diffusion coefficient weighted by transport theory forward and adjoint flux solutions for the heterogeneous cell; and they are accompanied by lattice-cell flux discontinuity factors at the interfaces between adjacent cells in the fuel assembly. The final fuel-assembly-homogenized whole-reactor diffusion equation has cross sections and diffusion coefficients that are given for each assembly by the assembly average of the cell-homogenized cross sections and diffusion coefficient weighted by the cell-homogenized diffusion theory forward and adjoint flux solutions for the cell-homogenized heterogeneous fuel assembly, and they are supplemented by fuel-assembly flux discontinuity factors at the interfaces between adjacent fuel assemblies. After the equations at each of the three levels – the forward and adjoint transport equations for each different type cell, the cell-homogenized forward and adjoint diffusion equations for each type assembly, and the assembly-homogenized diffusion equation for the whole reactor – are solved numerically, these solutions can be used to reconstruct the spatially detailed leadingorder transport theory solution for the original heterogeneous reactor within the theoretical framework provided by the results of the three-spatial-scale asymptotic expansion. The simple recipe for doing this follows directly from the asymptotic analysis. The numerical solutions for the assembly-averaged fluxes generated in the whole reactor nodal diffusion calculation are simply multiplied by the normalized numerical solutions for the cell-averaged fluxes produced in a nodal diffusion calculation for the corresponding assemblies, and then the resulting cell averages are multiplied by the normalized numerical solutions for the heterogeneous cells that result form the transport theory calculation for the corresponding cells. Thus, the results of the three-spatial-scale asymptotic expansion provide not only lattice-cell and fuel-assembly homogenized diffusion equations for coarse-mesh and nodal calculations, but also a self-consistent framework for the detailed reconstruction of the approximate transport theory flux at the spatial level of individual fuel pins [81]. Results of nodal diffusion criticality calculations, for rather realistic and detailed models of a pressurized water reactor, using fuel assembly size nodes, based on each of these three homogenization procedures have been reported and they were in good agreement with the results of the fine-mesh calculations for the corresponding original heterogeneous reactor models. The calculations carried out

438

J. Dorning

in the framework of the transport-theory-based two-scale expansion [80] were done using the discrete nodal transport method [70] for the two-dimensional heterogeneous fuel assembly calculations and the nodal Green’s function method [64] for the assembly-homogenized whole-reactor diffusion calculation; and those carried out in the framework of the diffusion-theory-based two-scale expansion [106] were done using the nodal Green’s method for both the two-dimensional heterogeneous fuel assembly diffusion calculations and the assembly-homogenized whole-reactor diffusion calculation. Finally, those carried out in the framework provided by the transport-theory-based three-scale expansion [81] were done using the discrete nodal transport method for the two-dimensional heterogeneous lattice cell calculations and the nodal Green’s function method for both the lattice-cell-homogenized fuel-assembly diffusion calculations and the fuel-assembly homogenized wholereactor diffusion calculation; and the resulting reconstructed pin powers were in very good agreement with those of the reference solution which was generated by a discrete nodal transport method calculation for the original heterogeneous wholereactor model [81]. These computations were done to validate the homogenization and flux reconstruction procedures that followed from the asymptotic expansion analyses and to explore their compatibility with coarse-mesh and nodal methods, for which they were developed. The calculations demonstrated that nodal diffusion methods are, indeed, very compatible with these homogenization procedures which, in fact, do provide a very suitable framework for their implementation. These computations, however, were done for a one-energy-group heterogeneous reactor model; and, although the extension of the asymptotic expansion development of homogenization procedures to two-group and few-group models is straightforward and has already been briefly reported [107], and a homogenization procedure for the multigroup equations has been developed based on a different approach [108], the implementation of two-group or few-group nodal diffusion methods in conjunction with the resulting homogenization procedures remains to be done. Further, the one-group calculations mentioned above were done in the context of reactor statics, i.e., criticality calculations, not reactor kinetics – which is the focus of this chapter. Thus, the implementation of these homogenization procedures – most importantly the transport-theory-based three-scale simultaneous lattice-cell and fuel-assembly homogenization procedure – has not yet been combined with one-group, two-group, or few-group practical reactor kinetics calculations. That remains for the future – hopefully, the near future. In this context, a few closing remarks on homogenization and reactor kinetics seem to be in order. Some years ago, a simultaneous development of lattice-cell homogenization and point kinetics was derived from the one-group time-dependent neutron transport equation and the coupled delayed-neutron precursor equations via a two-spatial-scale asymptotic expansion [109]. This expansion was not in the ratio of the prompt neutron lifetime to a characteristic precursor decay time like the derivation of point reactor kinetics as an asymptotic approximation to timedependent transport theory mentioned in Section 8.3.3 above. Rather, it was in the same small parameter used to develop the transport-theory-based homogenization

8


439

procedures just described – the ratio of the scattering mean free path to the characteristic reactor dimension – which, of course, is essential if the derivation is to lead to homogenized cross sections and diffusion coefficients. The end result of that asymptotic development was point kinetics equations based on forward and adjoint shape functions – the forward and adjoint fundamental eigenfuctions – obtained from instantaneous critical lattice-cell homogenized diffusion equations with cross sections and diffusion coefficients homogenized over lattice cells using forward and adjoint transport theory solutions for the heterogeneous cells. This derivation led to the point kinetics equations rather than the time-dependent diffusion equation because a slow time scale was used. If a less slow time scale had been introduced, a timedependent lattice-cell homogenized diffusion equation would have been the result. Clearly, this asymptotic development could be extended – by introducing three spatial scales – to arrive at a self-consistent derivation of lattice-cell and fuel-assembly homogenization procedures and point kinetics equations for use in point reactor kinetics calculations, and – using a less slow time scaling – to arrive at analogous homogenization procedures and cell- and assembly-homogenized whole-core diffusion equations for use in space–time reactor kinetics calculations. Further, if some of the ideas, to which reference was made above, were used it would be possible to start the two developments from a two-group or few-group transport theory description of a fully heterogeneous LWR and self-consistently derive both two-group or few-group time-dependent cell- and assembly-homogenized diffusion equations, and point reactor kinetics equations with cell- and assembly-homogenized forward and adjoint shape functions given by the solution to the associated whole-reactor few-group or two-group cell- and assembly-homogenized instantaneously critical diffusion equations. These equations, of course, would result, respectively from introducing a slow time scaling and a very slow time scaling. The final set of the three self-consistent models then would be: 1. The original few-group or two-group time-dependent transport equations for the fully heterogeneous reactor 2. The derived few-group or two-group cell- and assembly-homogenized timedependent diffusion equations for the whole reactor 3. The derived point kinetics equations for the reactor And the three self-consistent models could be used in a single reactor kinetics calculation; the first during time periods in the transient during which the time scale is fast, i.e., when the time derivative of the neutron flux is large; the second during periods when the time scale is intermediate, i.e., when the time derivative is neither large nor small; and the third, when the time scale is slow, i.e., when the time derivative is small. Moreover, since the asymptotic expansions, used to develop the second and third models, include correction terms – the first-order terms in the asymptotic approximations to the transport theory solution – these could be estimated and used to determine when, during a transient, a switch should be made from one model to another. Such an adaptive model calculation would be somewhat analogous to adaptive grid calculations which have long been used in computational fluid dynamics. There, time-step control is used to capture the correct dynamics

440

J. Dorning

when time-derivatives become large and spatial grid refinement is used to obtain the correct spatial variations when gradients become large – all within the context of a single model. Transitioning from one of the three reactor kinetics models described above to another of them in an “adaptive model calculation” also would result in the accurate calculation of various spatial gradients of the neutron flux, while the accurate calculation of transients associated with time derivatives of varying magnitude was being achieved. Interestingly, a calculation, not unrelated in general spirit, to the adaptive-model procedure just described, has recently been reported [82,83]. Although the computational models used in that work did not result from a single self-consistent derivation of homogenized diffusion equations and point reactor kinetics equations as approximations to the transport equation for a heterogeneous reactor, it nevertheless does represent the extension from variable-time-step adaptive-grid calculations to the important concept of adaptive model calculations. Thus, it will be summarized briefly in the next section.

8.5.5 Adaptive-Model Kinetics Calculations Adaptive grid methods have been used in conjunction with variable time-step methods in computational fluid dynamics for a few decades now. There, spatial grids are refined when necessary during flow simulations to accurately calculate the large gradients that develop in the vicinity of a shock-wave front in a compressible fluid flow calculation or in the vicinity of a weather front in a meteorology calculation. Variable model calculations also have been used in the context of compressible gas dynamics where, for example, the full nonlinear Boltzmann equation of the kinetic theory of gases is used in regions containing shock waves, boundary layers, etc., and hydrodynamics models derived as asymptotic approximations to the Boltzmann equation are used in the other – typically much larger – regions of the computational domain. The logical extension of adaptive grid methods to adaptive model methods may seem obvious; nevertheless, the nuclear reactor community waited a long time for its appearance. Happily, an adaptive model method has appeared – but only during the past decade, and just barely within the twentieth century (in 1999) [82, 83]. In that work, three models were used, but they were not the three models discussed here at the end of Section 8.5.4. Rather, they were a three-dimensional detailed diffusion theory model, a one-dimensional diffusion theory model and a point reactor kinetics model – which can be regarded as a “zero-dimensional” model. Both the one-dimensional diffusion theory model and the point kinetics model were derived from the three-dimensional diffusion theory model using shape functions based on the adiabatic approximation [31]; a three-dimensional shape function from the three-dimensional diffusion model was used in the point kinetics model, and a two-dimensional (transverse) shape function was used in the one-dimensional diffusion theory model. Since three-, one- and zero-dimensional models were used

8


441

this work was described as dimensionally adaptive reactor kinetics, rather than adaptive model reactor kinetics. The numerical method that was used to compute the space dependence of the neutron flux and precursor concentrations was based on the nodal expansion method – the first of the transverse-integration nodal methods developed for the diffusion equation (as mentioned in Section 8.5.3 above) [66, 67]. Hence, the numerical solutions to the equations of the three models were based on a unified approach. The criteria for changing, during a kinetics calculation, from a higherdimensional model to a lower-dimensional model were based on the change in time of the spatial flux distribution, and those for transitioning from a lower to a higher dimensional model were based on error estimates of the global neutron flux and of the reactivity. This dimensionally adaptive reactor kinetics method was combined with a widely used thermal hydraulics model and applied to three postulated reactor safety transients: the hot zero-power main-steam-line break, the rapid ejection of a single control assembly, and the boron dilution transient after a small-break lossof-coolant accident. The calculations led to good agreement between the numerical results of the dimensionally adaptive reactor kinetics simulations and those of the full three-dimensional simulations, whereas the straightforward one-dimensional diffusion theory simulations and the straightforward point kinetics simulations led to errors that were unacceptable. The results that follow from this study demonstrate the potential value of the adaptive model approaches to reactor kinetics and reactor dynamics calculations discussed at the end of Section 8.5.4 above.

8.6 Reactor Dynamics It is not uncommon to use the term “reactor kinetics” in reference to the time dependence of the neutron population in a reactor in which cross sections, etc. do not depend upon the neutron flux through temperature feedback or other feedback mechanisms – and the equations are therefore linear. Conversely, the words “reactor dynamics” are frequently employed in connection with reactor transients in which the cross sections, etc. do depend upon the flux via temperature, or other, feedback mechanisms – and, of course, the equations are nonlinear. It is precisely this usage that will be adopted here. As one might easily imagine, reactor dynamics is therefore a very broad subject, and it can also be very specific to a specific type reactor since different nonlinear feedback mechanisms may dominate over the others for different types of reactors – pressurized water thermal reactors (PWRs), boiling water thermal reactors (BWRs), liquid metal fast reactors (LMFRs), etc. Thus, a general systematic discussion of reactor dynamics like the general systematic discussion of reactor kinetics above is not appropriate – nor is it even possible. Therefore, just a few special topics will be discussed briefly in this short section on reactor dynamics. The feedback mechanism that has been of broadest historical interest in reactor dynamics is probably the Doppler effect, or more precisely the Doppler broadening of the giant narrow resonances in neutron cross sections which are well described

442

J. Dorning

by the Breit–Wigner single-level resonance formula in the resonance region through which neutrons slow down. This effective broadening of the narrow resonances results, as temperature increases, from the increased thermal motion of the atoms containing the nuclei with which the neutrons interact in the narrow resonances. The broader thermal distribution of the nuclei with increasing temperature leads to an increased number of neutrons having relative velocities that are in the narrow energy interval associated with the giant resonance than would be the case if the nuclei were at rest. Thus, the neutrons, which undergo slowing down via discrete energy decrements, are less likely to escape absorption in an absorption resonance. Hence the resonance absorption escape probability is decreased as a result of this Doppler broadening effect. (For a useful, but somewhat dated, discussion of the broadening of overlapping resonances and the unresolved resonance region, see [7]). Of course, there also are scattering and fission resonances; hence, the net effect of Doppler broadening could be either to increase or decrease the reactivity that results from the temperature rise associated with an increase in reactor power. When the result is an increase in the reactivity, the Doppler feedback is said to be positive; when it is a decrease, the Doppler feedback is said to be negative. In heterogeneous (real) reactors there are other important factors, such as coolant temperature, etc. that can also make important contributions to temperature feedback. (For some useful discussion of some of these factors, that is of historical interest and still relevant, see [13]). Of course, for obvious safety reasons, most power reactors are designed so that they have negative temperature feedback reactivity, making their control easier – although some Soviet-era thermal reactors had positive temperature feedback under certain conditions. With the advent of large light-water thermal reactors in the 1960s concern arose over the possibility of space-dependent xenon oscillations in these large systems. This worry arose because the fission product xenon-135, whose half-life is 9.2 h, has a very large thermal neutron absorption cross section of about three million barns; and because of this, in addition to the well-known poisoning effect of xenon135, there was a fear that large thermal reactors operating at high flux levels might undergo local oscillations in power as a result of local perturbations leading to local flux increases resulting in local xenon-135 burn-up, etc. [110–112]. During the late 1950s and throughout the 1960s, these problems and other problems in nonlinear reactor dynamics were the subjects of some very exciting research. One of the pioneers of this research was Brookhaven National Laboratory’s Jack Chernick – whose photo (from [113]) appears on the next page. He was “doing” nonlinear dynamical systems, phase portraits, “and all that” applied to nuclear reactor analysis long before nonlinear dynamics became so fashionable as it is today – and long, long before I knew anything about the subject. He also authored a lovely review on nonlinear reactor dynamics problems over 40 years ago [114]! This article also appears in the selection of his works, Jack Chernick 1911–1971: Selected Papers [113]. (Although spatial xenon oscillations in large thermal reactors turned out not to be a serious problem because, due to their slow time scale, reactor operators are able to control them, the pioneering research done on that problem and many other problems in nonlinear reactor dynamics during that era has provided the foundations

8


443

The late Jack Chernick. (With permission from Brookhaven National Laboratory.)

for much of the more recent application of bifurcation theory, modern nonlinear dynamical systems, and deterministic chaos to nonlinear problems that are of considerable importance to nuclear reactor dynamics – such as the stability of boiling water reactors (BWRs).) It seems appropriate that this brief discussion of a few selected topics in reactor dynamics end with a short summary of some aspects of the nonlinear dynamics of BWRs – since it is a very interesting application of bifurcation theory and nonlinear dynamical systems to reactor analysis, and the stability of BWRs is certainly an important problem. From a thermal hydraulics perspective, a BWR is basically an assemblage (a horizontal array) of vertical heated channels, each with subcooled fluid (water) entering at the bottom, reaching saturation temperature as it flows upward, and exiting the top as a two-phase liquid–vapor (steam–water) mixture. This flow, of course, is in parallel, and, because there are a large number of these parallel channels, each channel is subjected to an external pressure drop along its vertical length that is essentially constant in time. Interestingly, instabilities of flow in such “boiling channels” were a concern long before the advent of BWRs [115]. (Related problems date back to the early days of refrigeration engineering.) A textbook introductory discussion of instabilities in two-phase flow in parallel heated channels in the context of the thermal hydraulics of BWRs is available [116]. It has long been known that one type instability in such a heated channel is the “flow excursion” or “Ledinegg” instability [115, 116]. This instability results from a system parameter (such as the external pressure drop or the heat flux supplied to the channel) being passed through a saddle-node bifurcation value. (A discussion of this bifurcation,

444

J. Dorning

also known as the turning point bifurcation or the limit point bifurcation, along with a discussion of Hopf bifurcation – which also is very relevant to the instability of flow in heated channels and BWRs – is available in a very readable introductory textbook on nonlinear dynamics and chaos [117]; and more advanced discussion of these and many related topics in bifurcation theory and modern nonlinear dynamical systems analysis can be found in a very well-written higher-level textbook [118]). The engineers who designed the original BWRs in 1950s were well aware of the Ledinegg instability, and they designed the reactors so that they could be operated in a regime that was not close in parameter space to the related stability boundary (or bifurcation set). However, there was another flow instability which could not be completely avoided by BWR design. Avoiding it requires a BWR operating strategy as well. This instability, which was not unknown at the time, corresponds to a transition from the steady flow operation of a heated two-phase flow channel, as a control parameter is varied, to growing “density-wave” oscillations in the two-phase mixture exiting the top of the channel. A theoretical analysis carried out using the homogeneous equilibrium model for the two-phase flow showed that the transition, from the stable steady-state or the stable fixed point of the system to an unstable fixed point and stable oscillations or a stable limit cycle – the density-wave oscillations, was the result of a supercritical Hopf bifurcation [119]. Subsequent analysis, using the more complicated and more realistic drift flux model for the two-phase flow, led to the same result [120]. Since a simple parallel two-phase-flow heated channel undergoes a Hopf bifurcation it is not surprising that a BWR – which although a much more complicated engineering system is nevertheless, as indicated above, basically an assemblage of such channels – might also undergo a Hopf bifurcation. This is indeed the case and such a transition of an operating BWR – from steady 100% power to power oscillations of about 18% above and below that average power – occurred at the LaSalle-2 plant near Chicago, IL in 1988 following a valving error made by one of the plant engineers. Physically, of course, the oscillation of the length of the two-phase region in the BWR channels causes an oscillation in the scattering cross section – due to the great difference between the density of water and steam – and therefore in the diffusion coefficient and overall neutron leakage from the reactor. This, in turn, causes the neutron flux and fission rate to oscillate, which causes the power and the heat flux to the coolant channels to oscillate, etc. (Among the other transitions, from stable steady operation to power oscillations that have occurred, are those that have taken place at the Caorso BWR in Italy in 1984 and at the Laguna Verde BWR in 1995.) Theoretical studies of the nonlinear dynamics and stability of flow in two-phase flow heated channels, naturally, led to analogous studies of the nonlinear dynamics and stability of BWRs using simple models developed from the models used to study the heated channels. These simple models were generated by first combining the point reactor kinetics equations with the equations that described the heated channel through void-fraction- and Doppler-feedback in the reactivity. Subsequently, these models were extended to space–time kinetics by combining spatial modal expansions of the coupled neutron diffusion equation and coupled delayed neutron precursor equations with the heated channel equations. One such model,

8


445

which combined the space–time modal kinetics equations with a simplified homogeneous equilibrium model for the heated channel – along with models for the heat conduction in the fuel pin and clad and in the gap between the pin and the clad – led to a fairly serious model which although somewhat complicated still permitted the numerical calculation of the stability boundary – or bifurcation set – in parameter space, which when crossed led to supercritical Hopf bifurcation and a corresponding stable limit cycle – or power oscillation behavior [121]. Further, it was possible in that work to map the stability boundary from the parameter space to the powerflow diagram for the system – which corresponds to the power-flow map used by BWR operating engineers to avoid operation in the region, in the power-flow BWR operating regime, that would result in power oscillations. A projection of the parameter space with the (indistinguishable) stability boundaries or Hopf bifurcation sets – for the models based on the space–time modal kinetics equations and on the point reactor kinetics equations, but otherwise identical – and the BWR power-flow operating curve (or 100% rod line), mapped to that projection space, is shown in Fig. 8.11a; and the power-flow map with the operating curve (the natural circulation curve plus the 100% rod line) and the mapped stability boundaries appear in Fig. 8.11b [121]. The mapped stability boundaries or Hopf bifurcation sets overlap the upper left “corner” of the intersection of the steeply sloping natural circulation curve and the gently sloping hundred-percent rod line – indicating that operation of a BWR in this power-flow regime will result in power oscillations. This is precisely the region of the power-flow map that regulatory agencies require BWR operating engineers to avoid. To supplement the calculations of the stability boundaries, numerical simulations were done in that research and the decay ratios (DRs) were computed and they are tabulated in Fig. 8.11b, where it is clear that they are consistent with the close proximity of operating point “6” (on the 100% rod line) to the stability boundary. Although the two stability boundaries reported in that article were essentially indistinguishable (in Fig. 8.11a and b here), the numerical simulations reported there indicated that point kinetics was non-conservative vis-à-vis space–time kinetics, in that these two different representations of the neutron kinetics in the otherwise identical BWR models led to trajectories in phase space that converged to a stable operating point when point kinetics was used, but did not when space–time modal kinetics – which includes the important effect of the first spatial harmonic – was used. Further it was noted there that, as burn-up increases, the stability boundary moves closer to operating point “6” on the 100% rod line, making it less stable. For detailed studies of power oscillations in an operating BWR, simulations must be done using large-scale BWR codes that provide models of the neutron kinetics and of the system thermal hydraulics at a much more detailed level; however, the interpretation of such large-scale code simulations and the understanding of their engineering implications are greatly enhanced when they are based on the fundamentals of elementary bifurcation theory and modern nonlinear dynamical systems analysis – and the use of more manageable models of the type just discussed which capture the salient features of these BWR systems.

446

J. Dorning

DIMENSIONLESS INLET SUBCOOLING NUMBER

a

3 100% ROD LINE SB FOR MODEL KINETICS SB FOR POINT KINETICS

2.5

2 7 STABLE

1.5 6 5

1

4

0.5

3

UNSTABLE

1

2

0 0

b

5 30 10 15 20 25 DIMENSIONLESS EXTERNAL PRESSURE DROP

1 100% ROD LINE NATURAL CIRCULATION + TSB FOR MODAL KINETICS TSB FOR POINT KINETICS

0.9 0.8 NORMALIZED POWER

35

0.7

1 2 3

4 6

0.6 7

+

UNSTABLE

0.5

5 STABLE

+ +

0.4 0.3

+

0.2

+

0.1 0 0

+

0.1

Point 1 2 3 4 5 6 7

DR 0.42 0.48 0.52 0.78 0.94 0.97 0.90

+ 0.2

0.3

0.7 0.5 0.6 0.4 NORMALIZED FLOW

0.8

0.9

1

Fig. 8.11 (a) A projection of the parameter space with the (indistinguishable) stability boundaries for point kinetics and for modal kinetics and the 100% rod line mapped from the power flow map. (b) The power flow map with the stability boundaries mapped from the projection of the parameter space (with permission of the American Nuclear Society, Copyright October 1997 by the American Nuclear Society, La Grange, IL)

8


447

The nonlinear feedback mechanisms mentioned in this short selection of topics in reactor dynamics – Doppler feedback, xenon feedback, and void-fraction feedback – are but three of many nonlinear feedback phenomena that are ubiquitous in the very broad area of reactor dynamics, which, of course, also includes the reactor dynamics of severe and mild accidents – to mention just two of many important applications. Useful, but unfortunately not very recent, discussions of these and other feedback mechanisms can be found in [7, 13, 14].

8.7 Epilogue: 1934–1999 (and Prologue: 2000 – and Beyond) “A Century in Review – A Century Anew!” Selected developments in reactor kinetics in the twentieth century – or, more precisely, the final two-thirds of the twentieth century – have been reviewed. The earliest works cited were the 1934 and 1936 articles by Enrico Fermi and his collaborators at the Università di Roma on the time-dependent neutron diffusion equation published in the pre-dawn of the nuclear era, and the latest major reference cited was published in 1999 – the sunset of the twentieth century. During this interval of two-thirds of a century – from the genesis of nuclear reactor kinetics, and nuclear reactor physics in general, to the present state-of-the-art – the discipline has evolved from analytical solutions to the one-speed neutron diffusion equation for bare homogeneous simple reactor geometries, and analytical solutions to the point reactor kinetics equations, to numerical simulations of neutron kinetics, and reactor dynamics with thermal hydraulic feedback, in detailed complex heterogeneous reactor configurations, using large-scale supercomputers and advanced workstations. This dramatic progress in capability, of course, has been greatly facilitated by the enormous advances in the development of digital computers. Nevertheless, in retrospect it is still quite dramatic! The preceding six sections of this chapter comprise a review, albeit a limited review, of some of the more important developments of nuclear reactor kinetics during the twentieth century. But the title of this chapter is “Nuclear Reactor Kinetics: 1934–1999 and Beyond.” So this epilogue, reflecting back on the developments of the twentieth century, should also address the “and Beyond” part of the title and comprise a prologue to the twenty-first century. A prediction of what the important developments – or even the important needs – of reactor kinetics and reactor dynamics in the next century, however, is well beyond the range of my crystal ball. Thus, a limited speculation comprising some thoughts, and related caveats, on a few near-term concerns that might arise, is all that will be proffered. But first, I shall begin this closing section with an observation on where we are and where we might go fairly easily and fairly quickly in an important area of computational space–time reactor kinetics.

448

J. Dorning

8.7.1 Adaptive-Model Reactor Kinetics The general capability that now exists in the area of space–time reactor kinetics – based on modern nodal methods combined with recently developed homogenization procedures, as well as other advanced computational procedures – puts us in a good position for carrying out many important analyses of the kinetics and dynamics of currently planned LWR systems and many new reactor concepts as well. However, it certainly would be very desirable to be able to use a detailed space–time few-group transport-theory description of reactor kinetics during portions of a simulation when that is necessary, and to switch to a few-group (or collapsed group) diffusion-theory description when that is adequate, and even to a quasi-static kinetics calculation and a point kinetics calculation during the same simulation when they are adequate. This, of course, means using an adaptive model procedure of the general type discussed at the end of Section 8.5.4 above. A self-consistent framework for such an adaptive model procedure can be developed by systematically deriving few-group or multigroup spatially homogenized diffusion equations for reactor kinetics calculations from the few-group or multigroup space–time transport equations – and also deriving quasi-statics and point kinetics equations from the same few-group or multigroup space–time transport equations – all within the same general analytical development. This development can be carried out by expanding the solution to the few-group or multigroup transport equations in an asymptotic expansion in the ratio of the mean free path to the characteristic overall reactor dimension, and introducing various scalings for the time derivative (and fixed source when one is present) in this small parameter, as described briefly in Section 8.5.4 above and as done in several of the references cited there. The result will be reactor kinetics equations, of the three types just mentioned, that will yield as solutions the leading-order term in the asymptotic solution to the full few-group or multigroup transport equations. Which type of equations result – the original heterogeneous space–time transport equations, homogenized space–time diffusion equations, or quasi-statics or point kinetics equations – of course will depend upon the time scaling used; however, all will have been derived in the same general development based on the same asymptotic expansion of the solution to the transport equations. Depending on the amount of within-group neutron scattering and withingroup fission neutron regeneration, some of the groups of the original few or many energy groups in the fundamental transport theory description may be collapsed in the resulting homogenized space–time diffusion theory equations, and in the homogenized space-dependent diffusion equations for the shape functions for the resulting quasi-statics and point kinetics equations. Because the second terms in the asymptotic expansions provide corrections to the leading-order asymptotic solutions, they could be estimated and used to determine – in a completely consistent manner – when, during a transient calculation, the switch to a higher-level model (a more precise description) or to a lower-level model (a less precise description) should be made in order to retain some overall desired accuracy in the adaptive model calculation.

8


449

In applications in which the mean free path is not small compared to the overall characteristic reactor dimension, such as certain fast reactors, the ratio of the prompt neutron lifetime to a typical delayed neutron precursor decay time could be used as the small parameter instead of the ratio of the mean free path to the characteristic reactor dimension. (An approach of this type also was mentioned briefly in Section 8.5.4 above.) This, naturally, would not lead to homogenized space–time diffusion equations as a mid-level description, but rather to the original multigroup space–time transport theory description as the higher-level model and quasi-statics or point kinetics with multigroup – possibly with some groups collapsed – transport theory shape functions as a lower-level model. Of course, this would be appropriate since a diffusion theory description would not be a good representation of the neutron kinetics of such a system. Finally, this systematic development of such adaptive model frameworks based on asymptotic expansions could be supplemented by variational formulations based on various straightforward extensions of the Roussopoulos functional discussed in Section 8.5.2 above. Although the asymptotic formulations have the distinct advantage of leading deductively to the forms of the various lower-level models, they have, as a result, the disadvantage of being less “flexible” than variational methods. The latter, which are constructive rather than deductive, allow the analyst to choose the trial functions and thereby influence the outcome, i.e., the forms of the various lower-level models. This flexibility can be advantageous; however, it also can be dangerous. Hence, a combination of asymptotic developments of the lower-level models, followed by variational formulations – which allow the engineer’s intuition to enter through the trial functions and achieve “fine-tuned” lower-level models – may be an appropriate course.

8.7.2 Reactor Dynamics of Advanced Reactors As stated in the first paragraph of this epilogue/prologue, only a few – actually, only a couple – thoughts on near-term concerns that might arise in the area of reactor dynamics will be offered here. New generation reactor concepts have been driven, to a significant extent, by concerns for reactor safety – especially the natural convection cooling of the reactor core following emergency shutdown. A response to this concern in the context of BWR technology has been the introduction of the advanced boiling water reactor (ABWR) and the evolutionary simplified boiling water reactor (ESBWR), in both of which the natural convection cooling following shutdown is greatly enhanced due to the addition of unheated riser sections above the core coolant channels. The additional buoyancy due to the vapor in the steam–water mixture flowing in these risers causes this enhanced convection and the resulting improved safety following emergency shutdown; however, not surprisingly, it has been shown that the addition of these risers also has a destabilizing

450

J. Dorning

effect, increasing the risk of a transition from stable steady-state operation to power oscillations via a Hopf bifurcation as described in Section 8.6 above for conventional BWRs [122]. Many other advanced reactor concepts are based on pool-type designs in which the enhanced convective cooling of the core results from the fact that it is immersed in a large pool of liquid (e.g., liquid sodium or molten lead). This clearly is a very desirable safety feature; however, the flow patterns in such heated pools are wellknown to be very complicated, and, due to the strongly nonlinear nature of these flows multiple steady states can exist, and dramatic transitions among them can occur as systems parameters – such as the heat source from the shutdown core – vary. Even in a simple rectangular cavity heated from below (the so-called Bénard convection problem) and in an equally simple rectangular cavity with spatially uniform volumetric heating, the possible flow patterns and the static and dynamic bifurcations that lead to transitions among them – as the system parameters (the magnitude of the heat source and the height to width ratio of the cavity) are varied – are very complicated [123]. They can be steady, periodic, quasi-periodic, or even chaotic [124]. Thus careful attention must be paid not only to the nonlinear dynamics of the strongly nonlinear convective flows in these pool-type configurations, but also to the nonlinear dynamics of the coupled systems comprising the reactor core and the complex convective flows in these pool-type configurations.

8.7.3 Reactor Dynamics in the Twenty-First Century The concerns sketched in the previous section with respect to the stability of ABWRs and ESBWRs in connection with the riser sections in their designs, and with respect to the many complex flow states that could occur in pool-type reactors are but two examples of the many potential problems that could arise in complex – or even simple – advanced reactor systems that are likely to be pursued during the twenty-first century. There are many other examples; due to the fundamental strong nonlinearity of many of these systems, bifurcation phenomena and complex nonlinear dynamical behavior are ubiquitous. Thus, as the twenty-first century unfolds, in order to solve the important problems that arise, nuclear reactor dynamics will have to incorporate more and more analytic – and especially numerical – bifurcation analysis in conjunction with modern nonlinear dynamical systems analysis. Happily, the state-of-the-art in these two areas has advanced sufficiently far during the last quarter of the twentieth century that it is readily available now at the dawn of the twenty-first century for very practical application to real-world problems in engineering and applied science – including the advanced analysis of nuclear reactor dynamics. Acknowledgments I am very grateful to Ms. Alice Rice of Oak Ridge National Laboratory for her expert preparation of this manuscript. The writing was done at Le Carlina Lodge in Biarritz, France, and I should like to thank Mlle. Juliette Bégué for arranging our stay there, and I especially want to thank M. Jean Bernes, whose

8


451

kind and gentle heart, and quiet efficiency, kept everything running so smoothly at Le Carlina during the writing. And I enthusiastically add my gratitude to Mike and Paula Dudley for providing just the right amount of distraction at Le Carlina. Les vagues were also great! Last, but not least, I thank Helen, my wife of 45 years, for her understanding and encouragement, not only during the writing of this chapter, but throughout our life together.

References 1. Duderstadt JJ, Hamilton LJ (1976) Nuclear reactor analysis. Wiley, New York 2. Henry AF (1975) Nuclear-reactor analysis. MIT Press, Cambridge, MA 3. Glasstone S, Edlund MC (1952) The elements of nuclear reactor theory. Van Nostrand, New York 4. Murray RL (1957) Nuclear reactor physics. Prentice-Hall, Englewood Cliffs, NJ 5. Meghreblian RV, Holmes DK (1960) Reactor analysis. McGraw-Hill, New York 6. Lamarsh JR (1966) Introduction to nuclear reactor theory. Addison-Wesley Series in Nuclear Engineering. Addison-Wesley, Reading, MA 7. Bell GI, Glasstone S (1970) Nuclear reactor theory. Van Nostrand Reinhold, New York 8. Weinberg AM, Wigner EP (1958) The physical theory of neutron chain reactors. The University of Chicago Press, Chicago, IL 9. Davison B (1958) Neutron transport theory. The International Series of Monographs on Physics. The Oxford University Press, London/New York 10. Case KM, Zweifel PF (1967) Linear transport theory. Addison-Wesley, Reading, MA 11. Duderstadt JJ, Martin WR (1979) Transport theory. Wiley, New York 12. Keepin GR (1965) Physics of nuclear kinetics. Addison-Wesley Series in Nuclear Science and Engineering. Addison-Wesley, Reading, MA 13. Hetrick DL (1971) Dynamics of nuclear reactors. The University of Chicago Press, Chicago, IL 14. Akcasu Z, Lellouche GS, Shotkin LM (1971) Mathematical methods in nuclear reactor dynamics. Academic, New York 15. Stacey WM (1969) Space-time nuclear reactor kinetics. Nuclear science and technology. Academic, New York 16. Segrè E (Chairman of the Editorial Board) (1962) Enrico Fermi – Collected Papers, vol I, Italy 1921–1938. The University of Chicago Press, Chicago, IL 17. Segrè E (Chairman of the Editorial Board) (1962) Enrico Fermi – Collected Papers, vol II, United States 1939–1954. The University of Chicago Press, Chicago, IL 18. Wightman AS (Ed) (1992) The collected works of Eugene Paul Wigner, Part A The Scientific Papers, vol V. In: Weinberg AM (ed) Nuclear energy. Springer, Berlin/New York 19. Amaldi E, Fermi E (1936) Sull’Assorbimento dei Neutroni Lenti. – III. Ricerca Scientifica 7(1):56–59 20. FermiE, Amaldi E, Pontecorvo B, Rasetti F, Segrè E (1934) Influence of hydrogeneous substances on the radioactivity produced by neutrons – I. Ricerca Scientifica 5(2):282–283 21. Fermi E (1936) Sul Moto dei Neutroni nelle Sostanze Idrogenate. Ricerca Scientifica. 7(2):13–52 22. Wigner EP (1942) On Variations of the power output in a running pile, CP-351 (Chicago Project Report), November 11, Chicago University, Metallurgical Laboratory: Chicago, IL (see in [18] above) 23. Fermi E (1942) Problem of time dependence of the reaction rate: effect of delayed neutrons emission, CP-291 (Chicago Project Report), October 7, Chicago University, Metallurgical Laboratory: Chicago, IL (see [17] above) 24. Parks DE, Nelkin MS, Wikner NF, Beyster JR (1970) Slow Neutron Scattering and Thermalization, with Reactor Applications. WA Benjamin, New York

452

J. Dorning

25. Beckurts KH, Wirtz K (1964) Neutron physics. Springer, New York 26. Williams MMR (1966) The slowing down and thermalization of neutrons. North-Holland, Amsterdam 27. Ussachoff LN (1955) Equation for the importance of neutrons, kinetics and the theory of perturbations. Proceedings international conference on peaceful uses of atomic energy, Geneva, 1955, P/656, 5: p. 503–510. Columbia University Press, New York 28. Henry AF (1955) Computation of parameters appearing in the reactor kinetics equations, WAPD-124, December 1955, Westinghouse Electric Corp. Bettis Plant, Pittsburgh, PA 29. Henry AF (1958) Application of reactor kinetics to the analysis of experiments. Nucl Sci Eng 3:52–70 30. Lewins J (1965) Importance, the adjoint function: the physical basis of variational and perturbation theory in transport and diffusion problems. Pergamon, Oxford/New York 31. Henry AF, Curlee NJ (1958) Verification of a method of treating neutron space-time problems. Nucl Sci Eng 4:727 32. Kaplan S (1961) The property of finality and the analysis of problems in reactor space-time kinetics by various modal expansions. Nucl Sci Eng 9:357 33. Henry AF (1964) The application of inhour modes to the description of nonseparable reactor transients. Nucl Sci Eng 20:338–351 34. Ott K, Meneley D (1969) Accuracy of the quasistatic treatment of spatial reactor kinetics. Nucl Sci Eng 36(3):402 35. Ott K, Madell JT (1966) Quasistatic treatment of spatial phenomena in reactor dynamics. Nucl Sci Eng 26(4):563 36. Stacey Jr, WM, Adams CH (1967) Time-integrated method: a quasistatic neutron space-time approximation. Trans Am Nucl Soc 10:251–252 37. Lewins J (1960) The approximate separation of kinetics problems into time and space functions by a variational principle. J Nucl Energy Part A 12:108 38. Gyftopoulos EP (1964) General reactor dynamics. In: Thompson TJ, Beckerly JG (eds) The technology of nuclear reactor safety, vol I. MIT Press, Cambridge, MA, pp 175–204 39. Becker M (1968) A generalized formulation of point nuclear reactor kinetics equations. Nucl Sci Eng 31(3):458–463 40. Dorning J, Spiga G (1978) Point kinetics as an asymptotic representation of spacedependent, energy-dependent, and angle-dependent reactor kinetics. Trans Am Nucl Soc 28(June):761–762 41. Stakgold I (1967) Boundary value problems of mathematical physics, vol 1. Macmillan, New York/Collier-Macmillan Ltd, London 42. Dorning JJ (1980) Point kinetics with omega-mode shape functions via multiple-time-scale asymptotics. Trans Am Nucl Soc 34(June):281–283 43. Hille E, Phillips RS (1957) Functional Analysis and Semi-Groups, Rev. Ed., American Mathematical Society, Providence, RI 44. Lehner J, Wing GM (1955) On the spectrum of an unsymmetric operator arising in the transport theory of neutrons. Commun Pure Appl Math 8:217–234 (see also Wing GM (1962) An introduction to transport theory. Wiley, New York) 45. Corngold N (1964) Some transient phenomena in thermalization. Part 1, Theory. Nucl Sci Eng 19:80 46. Bednarz R (1965) Spectrum of the Boltzmann operator with an isotropic thermalization kernel. In: Proceedings IAEA symposium on pulsed neutron research, vol I. IAEA, Vienna, p 259 47. Albertoni S, Montagnini B (1965) Some spectral properties of the transport equation and their relevance to pulsed neutron experiments. In: Proceedings IAEA symposium on pulsed neutron research, vol I. IAEA, Vienna, p 239 48. Dorning JJ, Thurber JK ( 1968) Unstable pulsed-neutron decay states in non-crystalline moderators. Trans Am Nucl Soc 11:290 49. Conn R, Corngold N (1969) A theory of pulsed neutron experiments in polycrystalline media. Nucl Sci Eng 37(1):85

8


453

50. Dorning JJ (1968) Time-dependent neutron thermalization in finite light water systems. Nucl Sci Eng 33(1):81–92 51. Nelkin M (1960) Scattering of slow neutrons by water. Phys Rev 119(2):741–746 52. Koppel JU, Young JA (1964) Neutron scattering by water taking into account the anisotropy of the molecular vibrations. Nucl Sci Eng 19:412 53. Küchle MF, Mitzel F, Wattecamps E, Werle H (1966) Proceedings international conference fast critical experiments and their analysis, ANL-7320. Argonne National Laboratory, Argonne, IL, p 506 54. Dorning J (1972) Fast and thermal prompt kinetics via models. In: Hetrick DL (ed) Dynamics of nuclear systems. The University of Arizona Press, Tucson, AZ 55. Dorning J, Nicolaenko B, Thurber JK (1969) Unstable decay states of far subcritical fast assemblies. Trans Am Nucl Soc 12(1):251 (see also Dorning J, Nicolaenko B, Thurber JK (1969) Some mathematical considerations on Cauchy integrals related to pulsed neutron dispersion relations. In: Neutron transport theory conference. ORO-3858–1, U.S. Atomic Energy Commission, Washington, DC, pp 76–121) 56. Dorning J, Nicolaenko B, Thurber JK (1969) Theory of pulse decay in subcritical fast systems with general cross sections. Trans Am Nucl Soc 12(2):656 57. Dorning J, Nicolaenko B, Thurber JK (1969) Time-behavior of pseudo-modes in subcritical fast systems described by exact elastic slowing down. Trans Am Nucl Soc 12(2):657 58. Mitzel F, Schroeter E (1969) Gepulste Experimente an der Schnellen Unterkritischen Anlage Karlsruhe (SUAK) und Vergleich der Ergebnisse mit Berechnungen. Nukleonik 12:110–117 59. Gelbard EM (1973) Twenty years of reactor computation in retrospect. In: Proceedings conference numerical models and computational techniques for analysis of nuclear systems. NTIS, US Department of Commerce, Springfield, VA, p xxi–xxvi 60. Burns TJ, Dorning JJ (1973) A partial current balance method for space-, and energydependent reactor calculations.” In: Proceedings national topical meeting on mathematical models and computational techniques for analysis of nuclear systems, vol II. CONF-73044. American Nuclear Society, LaGrange Park, IL, p VII.162–VII.178 (see also Burns TJ, Dorning JJ (1975) Multidimensional applications of an integral balance technique for neutron diffusion computations. In: National topical meeting on computational methods in nuclear engineering, vol II. CONF-750413. US Energy Research and Development Agency, Washington, DC, pp V.57–V.68 61. Shober, RA, Henry AF (1976) Approximate analytical method for determining nodal fluxes. Trans Am Nucl Soc. 24(Nov 19):192–193 62. Birkhofer A, Werner W (1973) Efficiency of various methods for the analysis of space time kinetics. In: Proceedings national topical meeting on mathematical models and computational techniques for analysis of nuclear systems, vol II. CONF-730414. American Nuclear Society, LaGrange Park, IL, p IX–31 63. Lawrence RD, Dorning JJ (1978) Nodal Green-function method for multidimensional neutron diffusion calculations. Trans Am Nucl Soc 28(June):248–249 64. Lawrence RD, Dorning JJ (1980) A nodal Green-function method for multidimensional neutron diffusion calculations. Nucl Sci Eng 76(2):218–231 65. Greenman G, Smith KS, Henry AF (1979) Recent advances in an analytic nodal model for static and transient reactor analysis. In: Proceedings national topical meeting on computational methods in nuclear engineering, CONF-790402. American Nuclear Society, LaGrange Park, IL, p 3.49 66. Finnemann H (1975) A consistent nodal method for the analysis of space-time effects in large LWRs. In: Joint NEACRP/CSNI specialist’s meeting on new developments in three dimensional neutron kinetics and review of kinetics benchmark calculations. Laboratorium für Reaktorregelung und Anlagensicherung, Garching (Munich), Germany 67. Finnemann H, Bennewitz F, Wagner MR (1977) Interface current techniques for multidimensional reactor calculations. Atomkernenergie 30(2):123–128 68. Wagner MR (1979) A nodal discrete-ordinate method for the numerical solution of the multidimensional transport equation. In: Proceedings American nuclear society topical meeting on

454

69.

70.

71.

72.

73.

74. 75. 76. 77. 78. 79.

80.

81.

82.

83.

84. 85. 86.

J. Dorning computational methods in nuclear engineering, vol 2. American Nuclear Society, LaGrange Park, IL, p IV–177 Dorning, JJ (1979) Modern coarse mesh methods, a development of the 70’s. In: Proceedings American nuclear society topical meeting on computational methods in nuclear engineering, vol 1. American Nuclear Society, LaGrange Park, IL, p III–1 Lawrence RD, Dorning JJ (1979) New course-mesh diffusion and transport methods for the efficient numerical calculation of multi-dimensional reactor power distributions. In: OECD/NEA-CRP specialists’ meeting on calculations of 3-dimensional rating distributions in operating reactors. OECD, Paris, p 383 (see also Lawrence RD, Dorning JJ (1980) A nodal integral transport theory method for multidimensional reactor physics and shielding calculations. In: Proceedings of the American nuclear society international topical meeting on 1980 advances in reactor physics and shielding. American Nuclear Society, LaGrange Park, IL, pp 240–251 Walters WF, Odell RD (1981) Nodal methods for discrete-ordinates transport problems in (X, Y) geometry. In: Proceedings American nuclear society topical meeting on advances in mathematical methods for the solution of nuclear engineering problems, vol 1. American Nuclear Society, Munich, Germany, p 115 Badruzzaman A, Xie Z, Dorning JJ, Ullo JJ (1984) A discrete nodal transport method for three-dimensional reactor physics and shielding calculations. In: Proceedings American nuclear society topical meeting on reactor physics and shielding, vol 1. American Nuclear Society, LaGrange Park, IL, p 170 Dorning JJ (1985) Nodal transport methods after five years. In: Proceedings of the American nuclear society international topical meeting on advances in nuclear engineering computational methods. American Nuclear Society, LaGrange Park, IL, p 412 Smith KS (1986) Assembly homogenization techniques for light water-reactor analysis. Prog Nucl Energy 17(3):303–335 Lawrence RD (1986) Progress in nodal methods for the solution of the neutron diffusion and transport equations. Prog Nucl Energy 17(3):271–301 Koebke K (1978) A new approach to homogenization and group condensation. In: IAEA technical committee meeting on homogenization methods in reactor physics. IAEA, Vienna LarsenEW (1975) Neutron-transport and diffusion in inhomogeneous-media.1. J Math Phys 16(7):1421–1427 Larsen EW (1976) Neutron-transport and diffusion in inhomogeneous-media.2. Nucl Sci Eng 60(4):357–368 Chiang RT, Dorning JJ (1980) A homogenization theory for lattices with burnup and nonuniform loadings. In: Proceedings of the American nuclear society international topical meeting on 1980 advances in reactor physics and shielding. American Nuclear Society, LaGrange Park, IL, p 319–329 Zhang H, Rizwan-uddin, Dorning JJ (1997) Systematic homogenization and self-consistent flux and pin power reconstruction for nodal diffusion methods: Part II. Transport Theory Stat Phys 26(4):433–468 Zhang H, Rizwan-uddin, Dorning JJ (1997) A multiple-scales systematic theory for the simultaneous homogenization of lattice cells and fuel assemblies. Transport Theory Stat Phys 26(7):765–811 Jackson CJ, Cacuci DG, Finnemann HB (1999) Dimensionally adaptive neutron kinetics for multidimensional reactor safety transients - I: new features of RELAP5/PANBOX. Nucl Sci Eng 131(2):143–163 Jackson CJ, Cacuci DG, Finnemann HB (1999) Dimensionally adaptive neutron kinetics for multidimensional reactor safety transients – II: dimensionally adaptive switching algorithms. Nucl Sci Eng 131(2):164–186 Wachspress EL (1966) Iterative solution of elliptical systems. Prentice Hall, Englewood Cliffs, NJ Clark M, Hansen KF (1964) Numerical methods of reactor analysis. Academic, New York Varga RS (1962) Matrix iterative analysis. Prentice-Hall, Englewood Cliffs, NJ

8


455

87. Greenspan H, Kelber CN, Okrent D (1968) Computing methods in reactor physics. Gordon & Breach, New York 88. Cadwell WR, Henry AF, Vigilotti AJ (1964) WIGLE: A program for the solution of the two-group space-time diffusion equations in slab geometry, WAPD-TM-416. Pittsburgh, PA, Westinghouse Bettis Atomic Power Laboratory. See also: Henry AF, Vota AV (1965) WIGL2: A program for the solution of the one-dimensional, two-group, space-time diffusion equations accounting for temperature, xenon, and control feedback, WAPD-TM-532, Westinghouse Atomic Power Laboratory, Pittsburgh, PA 89. Yasinsky JB, Natelson M, Hageman LA (1968) TWIGL – a program to solve the twodimensional, two-group, space-time neutron diffusion equations with temperature feedback, WAPD-TM-748, Westinghouse Bettis Atomic Power Laboratory, Pittsburgh, PA 90. Courant R, Hilbert D (1953) Methods of mathematical physics. Interscience Publishers, New York 91. Kantorovich LV, Krylov VI (1964) Approximate methods of higher analysis. Interscience Publishers, New York 92. Kahan T, Rideau G (1952) Sur la Déduction de Divers Principes Variationnelles de la Théorie des Collisions a` Partir d’un Principe Unique. Journal de Physique et Radiation 13:326–332 93. Roussopoulos P (1953) Méthodes Variationnelles en Théorie des Collisions. Comptes Rendus Académie Science 236:1858–1856 94. Selengut DS (1959) Variational analysis of a multi-dimensional system. In: Hanford Works Quarterly Physics Report HW-59126. Hanford, WA, p 89–124 95. Stacey WM (1967) Variational functionals for space-time neutronics. Nucl Sci Eng 30(3):448 96. Michael EPE, Dorning JJ (2001) Studies on nodal methods for the time-dependent convection diffusion equation. In: Proceedings of the American nuclear society international meeting on mathematical methods for nuclear applications. American Nuclear Society, LaGrange Park, IL 97. Michael EPE, Dorning JJ (2001) A primitive-variable nodal method for the time-dependent Navier–Stokes equations. In: Proceedings of the American nuclear society international meeting on mathematical methods for nuclear applications. American Nuclear Society, LaGrange Park, IL 98. Mays BE, Dorning JJ (2001) A nodal integral method for the dissipative shallow-water equations with Coriolis force. In: Proceedings of the American nuclear society international meeting on mathematical methods for nuclear applications. American Nuclear Society, LaGrange Park, IL 99. Azmy YY, Dorning JJ (1983). A nodal integral approach to the numerical solution of partial differential equations. In: Proceedings American Nuclear Society International Topical Meeting on Advances in Reactor Computations, vol II. American Nuclear Society, LaGrange Park, IL, pp 893–909 100. Ougouag AM, Dorning JJ (1980) A hybrid nodal Green’s function method for multigroup multidimensional neutron diffusion calculations. In: Proceedings American nuclear society international topical meeting on 1980 advances in reactor physics and shielding. American Nuclear Society, LaGrange Park, IL, pp 228–239 101. Ougouag AM, Rajic HL (1988) Illico-HO: a self-consistent higher-order coarse-mesh nodal method. Nucl Sci Eng 100(3):332–341 102. Azmy YY (1988) The weighted diamond-difference form of nodal transport methods. Nucl Sci Eng 98(1):29–40 103. Koebke K (1981) Advances in homogenization and dehomogenization. In: Proceedings international topical meeting on advances in mathematical methods for the solution of nuclear engineering problems, vol 2. American Nuclear Society, Munich, Germany, p 59 104. Smith KS, Henry AF, Lorentz R (1980) Nucl Sci Eng 17:303–335 105. Grad H (1958) Principles of the kinetic theory of gasses. In: Flügge S (ed) Handbuch Der Physik, Band XII. Springer, Berlin 106. Zhang H, Rizwan-uddin, Dorning JJ (1995) Systematic homogenization and self-consistent flux and pin power reconstruction for nodal diffusion methods – 1: diffusion equation-based theory. Nucl Sci Eng 121(2):226–244

456

J. Dorning

107. Dorning J (2003) Homogenized multigroup and energy-dependent diffusion equations as asymptotic approximations to the Boltzmann equation. Trans Am Nucl Soc 89:313–316 108. Larsen EW (1977) A homogenized multigroup diffusion theory for the neutron transport equation. In: Glowinski R, Leon JL (eds) Third international symposium on computing methods in applied sciences and engineering. Lecture Noted in Mathematics, vol 704, 1979. Springer, Berlin, p 357 109. Dorning JJ, Chiang RT (1980) A simultaneous development of point kinetics and lattice homogenization theory. In: Proceedings American Nuclear Society International Topical Meeting on 1980 Advances in Reactor Physics and Shielding. American Nuclear Society, LaGrange Park, IL, pp 309–318 110. Chernick J (1960) The dynamics of a xenon-controlled reactor. Nucl Sci Eng 8:233–243 111. Chernick J, Lellouche G, Wollman W (1961) The effect of temperature on xenon instability. Nucl Sci Eng 10:120–131 112. Canosa J, Brooks H (1966) Xenon-induced oscillations. Nucl Sci Eng 26:237 113. Christoffersen AP, Levine MM, Michael PA, Price GA (1973) Jack Chernick 1911–1971– Selected Papers. Brookhaven National Laboratory, Upton, New York 114. Ledinegg M (1938) In stability of flow during natural and forced convection. Die Wärme 61(2):8 115. Chernick J (1962) A review of nonlinear reactor dynamics problems, BNL 774 (T-291). Physics: Reactor Technology – TID-4500, 18th edn. Brookhaven National Laboratory, Upton, New York 116. Lahey Jr, RT, Moody FJ (1977) The thermal-hydraulics of a boiling water nuclear reactor. The American Nuclear Society, LaGrange Park, IL 117. Strogatz SH (1994) Nonlinear dynamics and chaos. Perseus Books (Westview Press), Cambridge, MA 118. Guckenheimer J, Holmes P (1983) Nonlinear oscillations, dynamical systems, and bifurcations of vector fields. Springer, New York 119. Achard J-L, Drew DA, Lahey Jr RT (1985) The analysis of nonlinear density – wave oscillations in boiling channels. J Fluid Mech 98:616 120. Rizwan-uddin, Dorning JJ (1986) Some nonlinear dynamics of a heated channel. Nucl Eng Design 93:1–14 121. Karve AA, Dorning JJ, Rizwan-uddin (1997) Out-of-phase power oscillations in boiling water reactors. In: Proceedings of the joint international conference on mathematical methods and supercomputing for nuclear applications, vol 2. American Nuclear Society, LaGrange Park, IL, pp 1633–1647 122. Rizwan-uddin, Dorning JJ (1992) Stability of advanced boiling water reactors. In: American Nuclear Society Proceedings of the 1992 National Heat Transfer Conference, HTC vol 6. American Nuclear Society, LaGrange Park, IL, pp 161–170 123. Decker WJ, Dorning J (1996) Bifurcations and unfoldings in natural convection. Trans Am Nucl Soc 74:163–165 124. Decker WJ, Dorning J (1996) Transitions to chaos in natural convection. Trans Am Nucl Soc 74:165–166

8


457

Professor John J. Dorning graduated from the US Merchant Marine Academy at Kings Point, New York in 1959 with a BS in Marine Engineering, after which he served as an engineering officer aboard a destroyer tender in the US Navy and then as a ship’s engineer on various US merchant vessels. He earned his doctorate in Nuclear Science and Engineering from Columbia University in 1967 where his research on theoretical and computational studies of pulsed neutron die-away and neutron thermalization in moderating media led to his being selected to receive the American Nuclear Society’s Mark Mills Award. After 2 years as a member of the Reactor Theory Group at Brookhaven National Laboratory, in 1969 he became group leader of the Reactor Physics Analysis Group while simultaneously continuing his reactor theory research into several aspects of neutron transport theory. A year later he moved to the University of Illinois where he remained on the faculty of the Nuclear Engineering Program for fourteen years. Then, in 1984 he became Whitney Stone Professor of Nuclear Engineering, of Engineering Physics, and of Applied Mathematics at the University of Virginia. Over the years his research interests have expanded to many other areas of mathematical physics and computational science, including reactor kinetics, reactor physics, particle transport theory, computational methods development, plasma physics, fluid dynamics, heat transfer, and nonlinear dynamical systems, and deterministic chaos. His dedication to teaching and research has led to a number of honors and awards. A partial list of these includes: elected Fellow of the American Nuclear Society in 1978; of the American Physical Society in 1982; of the American Association for the Advancement of Science in 1986; awarded the American Society for Engineering Education’s Glenn Murphy Award for education in 1988; awarded the American Nuclear Society’s Arthur Holly Compton Award in 1998 recognizing his teaching and research accomplishments; awarded the ANS’s Eugene P. Wigner Reactor Physicist Award in 1999 recognizing his research contributions to reactor physics; and awarded the ANS’s Glen T. Seaborg Medal in 2002 recognizing his research contributions to the development of the peaceful uses of nuclear energy. In 1990 he received the US Department of Energy’s prestigious Ernest O. Lawrence Memorial Award for research; and in 2007 he was elected a member of the US National Academy of Engineering. He has been the editor of Transport Theory and Statistical Physics since 1995. He has authored and coauthored over 250 research articles on a wide variety of topics; and he has had the good fortune to have advised over 30 Ph.D. students, five of whom were selected to receive the prestigious ANS’s Mark Mills Award, and four of whom were honored by the University of Virginia’s coveted Allen Talbert Gwathmey Memorial Award.

Index

A Absolutely constant errors, 295 Acceleration methods diffuse-synthetic acceleration, 2, 17, 18, 20, 49, 54, 57–59, 71, 72 fission source acceleration, 2, 62 linear multifrequency-grey, 2, 17, 18, 60, 62 neutron thermal upscatter source acceleration, 62 upscatter acceleration, 2 Adams, M.L., 54 Adaptive Monte Carlo, 154. See also Sequential Monte Carlo Adjoint functions, 316–318, 323, 328, 329, 335, 336, 338, 342, 344–346, 349 Adjoint operator, 316, 330 Adjoint sensitivity analysis procedure global, 349–350 local, 291, 316–318, 321–346, 349, 350 Adjoint sensitivity system, 328, 336, 342, 345, 346, 349 Adkins, C.R., 232 Adler, D.B., 226–229 Adler, F.T., 226–229 Alcouffe, R.E., 59, 62 Algorithm, marching, 86 Amaldi, E., 377 Amster, H.J., 156 Analytical broadening line-shape functions evaluations Cauchy integral, 236, 237 continued fraction, 236 generalized, 234–235 traditional, 234 Andres, T.H., 128, 313

Angular flux, 5, 13, 20, 26, 34, 43, 50, 52, 55, 58, 59, 71, 72, 88, 89, 110 cell-average, 9–11, 24 cell-edge, 9, 36 Angular intensity, 15 Angular quadrature sets companion quadrature weights, 42 double Radau, 41 Galerkin quadrature, 2, 40–45, 48, 49, 75 Gauss–Legendre, 9 Lobatto, 41 Associated Legendre polynomials, 97 Asymptotic frequency, 135 ATTILA, 22 Azmy, Y.Y., 62, 375 B Balance equation, 9, 13, 15, 23, 26 Becker, M., 394 Bell, G.I., 411 Bell, G.L., 266 Benchmark, Azmy, 103 Bergstrom, V., 125 Bettonvil, B., 313 Bhat, M.R., 237 Bickley, W.C., 265, 270 Bilinear concomitant, 323, 345 Bode, H.W., 314 Boltzmann equation, 86 Boltzmann–Fokker–Planck equation, 18, 19 Boltzmann’s constant, 15 Boolean transformation, 95, 105, 111 Boundary condition, 89–93, 99, 106 essential, 89, 90, 92, 93, 110 incoming flux, 90, 106 modified natural, 90–93 natural, 89, 106 reflected, 90, 92, 93, 106 vacuum, 90, 92, 93, 95, 106, 112 zero current, 93

459

460 Boundary layers, 17, 20, 35, 36 Box, G.E.P., 311 Bramblett, E.T., 240 Breit, G., 218, 219, 224, 225, 230, 232, 246–248, 250, 251, 262, 267, 272, 274, 276 Brute-force method, 314, 348 Bubnov–Galerkin method, 424, 427 Buffon, G.C., 117

C Cacuci, D.G., 317, 318, 320, 325, 339, 350 Caflisch, R.E., 131 CANDU, 175, 189, 207–208, 215 Carlson, B.G., 1 Central limit theorem, 134, 296 Chandrasekhar, S., 14 Charged-particle transport, 2, 6, 14, 18–21, 30, 35, 45, 49, 75 Chebychev, 48 Chelson, P., 135, 136, 149 Chernick, J., 218, 257, 442, 443 Chicago Project Report CP-351, 379 Clement, R.T., 304 Cohen, E.R., 248, 249 Collision density, 144, 146, 147, 151 Collision probability black limit approximation, 264 convex geometries, 259 escape probability, 263–265 general features, 258–260 Nordheim’s approximation, 264 rational approximation, 259–261, 263, 265–266 theory, 111, 112 Computational effort, 311, 313, 317, 346, 348 Conditionally constant errors, 295 Confidence interval, 134, 140 Conjugate gradients, 96 Conover, W.J., 305 Conservation, neutron, 87, 107 Continuity conditions, 105 Continuous slowing-down approach, 253–254 Control variates, 142–146 Convergence rate, 129, 131, 138, 139 Corngold, N.C., 244, 264, 271, 406–409 Correlated parameters, 301, 304, 306, 347, 348 Correlated sampling, 145, 148, 151 Correlation, 119, 123, 125, 127, 128, 137, 143, 152 Correlation ratio-based method, 302, 308 Cotter, S.C., 312 Coulibaly, I., 131

Index Covariance, 298–300, 319 Coveyou, R.R., 122 Criticality, 355–369 Criticality safety, 355–359 Critical point, 302, 317, 324–329, 331, 337, 339, 350 Cross sections, 3–7, 9, 15–20, 31, 33–35, 38, 43–45, 48, 49, 62, 71, 105, 112, 316, 331 differential scattering, 3, 18, 43 extended transport correction, 45 group average Legendre flux moment, 50 Legendre moments, 3, 20, 39, 41, 43, 44, 59 restricted momentum transfer, 18 restricted stopping power, 18, 50 restricted transport-corrected scattering, 18 Rosseland-averaged, 17 scattering, 87, 96, 97 within-group, 6, 17, 60 Cruz, J.B., 315 Cukier, R.I., 308 Cullen, D.E., 283 Current, 86–88, 93, 110 Curvilinear geometries, 7 angular derivative, 2, 9, 36–37 conservative form, 8 flux dip, 36–37 starting direction, 9, 36–37

D Dancoff, S.M., 263, 264 Davison, B., 169, 171 Deif, A.S., 315 Derivative Gâteaux, 321, 322, 327, 333, 337, 339 partial, 314, 321, 327, 333, 337, 339 deSaussure, G., 228–230 Deterministic sensitivity analysis, 291, 302, 314–330, 347 Differential equations, 314, 315, 343–346 Differential, Gâteaux, 317, 318, 320, 348, 349 Diffusion particle, 331–339 theory, 85, 86, 92–96, 103, 110 Diffusive, 10, 11, 13, 18, 44, 51, 74 asymptotic diffusion limit, 2, 35, 36 optically-thick diffusion limit, 14, 31 solvability condition, 33 thick diffusion limit, 17, 30–36 Dirac functional, 324 Direct effect term, 321, 333, 344 Direct measurements, 297

Index Discontinuity factors, 192, 199, 206 Discrepancy, 125, 129, 135, 137, 138 Discrete ordinates, 1–75 approximation, 86, 96, 100–103, 107, 108 simplified, 103 Sn or SN , 1–75, 96, 102–104, 358–360, 362–363, 366, 368, 369, 371 Discrete-to-moment matrix, 40, 42 Discretization, 9, 18, 59, 72 characteristic, 2, 23–25 corner balance, 27 diamond-difference, 8, 10, 34, 51, 59 discontinuous finite element, 2, 20, 25, 54 finite difference, 1, 13, 14, 23 linear discontinuous, 25–27, 37, 43, 49–53, 71 bi-linear discontinuous, 34 lumped LD, 18, 27 trilinear discontinuous, 34, 51, 53 negative flux fixup, 10 nodal, 22, 27–29, 51–54 constant-constant nodal, 29 linear-linear nodal, 29 step characteristic, 34 weighted-diamond, 23, 25 Distribution, 118, 123, 125, 126, 134 joint probability, 298 moments, 297 normal, 296 Dolling, G., 233 Doppler-broadening crystalline binding effects, 232–233 practical kernels, 231–233 Dorning, J.J., 395 Douglas, C.H., 359 Dresner, L., 218, 286 Duality, 147, 148 Dual method, 87, 89–91, 112 Dunford, C.L., 240 Dunker, A.M., 317 Dyson, F.J., 274–276

E Eddington, A.S., 14 Egelstaff, P.A., 232 Eichenauer, J., 127 Eisenbud, L., 221–224 Elementary errors, 294, 295 Energy-balance equation, 15 Epistemic uncertainty, 293 Equations coupled, 87, 88, 97 differential, 314, 315, 343–346

461 differential–integral, 85 even parity, 86, 88, 90–92, 96, 98, 101, 102, 106, 109, 112 first order, 85–87, 96, 97, 99 odd parity, 89, 97 Ricatti, 339–342 second order, 85–88, 90, 91, 93, 96, 98, 99, 106, 112 sparse matrix, 112 Equidistribution, 125, 126 Ericson, T., 272 Error analysis, 118, 120, 129, 132–141 pseudorandom, 120, 129, 132–140 quasi-random, 129, 133–141 Error bound pseudorandom, 129, 140 quasi-random, 129, 134, 141 Error reduction, 120, 141–155 Errors absolutely constant, 295 conditionally constant, 295 elementary, 294, 295 purely random, 295, 296 quasi-random, 296 random, 294–296 systematic, 294, 296 Eslami, M., 315 Estimator, 135, 142, 148, 149, 152, 157. See also Random variable Euler–Lagrange equation, 92, 93, 99–101, 106, 111 Expected value, 118, 142, 143, 150–151, 156, 157. See also Mean Expected value estimator, 118, 142, 150, 157 Experiment data, 295, 296, 301, 349 one-at-a-time, 302, 311, 312

F Faure, H., 131 Fermi, E., 375–383, 447 Fiacco, A.V., 315 Finite difference, 85, 86 Finite element method, 86, 92, 94–96, 105, 106, 110–112 Finite elements, 201–203 dual formulation, 202 mixed formulation, 203 primal formulation, 202 First-order reliability method, 302, 307 Fischer, R.A., 310 Fission, 96

462 Flux angular incoming, 89 even parity, 101, 105, 108 odd parity, 87, 106, 108 of particles, 331–339 Fokker–Planck equation, 35, 49 Formal adjoint operator, 322, 334, 335, 341, 345 Forward sensitivity analysis procedure, 316–321, 329, 341, 344, 346, 348 Forward sensitivity model, 320 Forward variational model, 320 Fourier amplitude sensitivity test, 302, 308, 348 Fourier analysis, 18, 54–58, 61 iteration eigenvalue, 56 spectral radius, 57 Fourier coefficients, 322 Frank, P.M., 315 Fredholm alternative theorem, 396 Fröhner, F.H., 221, 237, 286 Functional, 297, 307, 314–317, 322–329, 332, 339, 343 Dirac, 324 even parity, 93 Heaviside, 336, 338 linear, 334, 341 nonlinear, 301, 317, 324 reduced form, 95, 100, 108 Fundamental lattice mode, 186 Fundamental mode, 253

G Gaseous diffusion plant, 356, 357, 359, 373 Gâteaux derivative, 321, 322, 327, 333, 337, 339 differential, 317, 318, 320, 348, 349 Gelbard, E.M., 196, 417 Generating sequences, 120, 122–134, 147 inversive congruential generator, 127–129 linear congruential generator, 123, 124, 126–128 linear recursive generator, 127 multiplicative congruential generator, 123 nonlinear congruential generator, 124, 127, 128 pseudorandom, congruential generator, 123, 124, 126–128 quasi-random, 124, 125, 129–131, 133, 141 Tausworthe generator, 128 Geometric convergence, 145 Ginsburg, M., 263, 264 Glasstone, S., 411

Index Global adjoint sensitivity analysis procedure, 349–350 Global analysis, 302, 349, 350 Goertzel, G., 254 Goldstein, R., 248, 249 Goodwin, E.T., 236 Goudsmit, S., 154 Gourevich, I.I., 218, 262 Greenman, G., 419 Greenspan, E., 349 Greuling, E., 254 Gyftopoulos, E.P., 394 H Hajas, W.C., 313 Halton, J.H., 130, 131, 136, 137, 140, 144, 145 Heat capacity, 15 Heaviside functional, 336, 338 Henrici, P., 236 Henry, A.F., 388, 389, 393, 394, 397, 419 Hickernell, F.J., 131, 146 Higher-order moments, 297, 300 Hilbert space, 321–323, 327, 334, 335, 341, 342 Hlawka, E., 130, 137–141 Holt, J., 212 Humblet, J., 230 Hummel, H.H., 253, 266 Hwang, R.N., 240 Hybrid finite element, 106 Hybrid method, 86, 106–112 Hybrid sequences, 131–133, 139–141 I Iman, R.L., 305 Implicit, 1, 4–7, 13, 14, 60 Importance function, 146–148, 157 Indirect effect term, 321, 323, 328, 333, 334, 336, 338, 341, 344, 345 Indirect measurements, 295–297, 349 Inner iteration, 18 Inner product, 321, 327–329, 334, 341 Integral method, even-parity, 86, 106–112 Interface conditions, 109 Iterative methods, 85, 96, 106, 110, 112 J Joint probability distribution, 298 K Kahn, H., 122 Kapur, P.L., 221, 223–225, 227–229

Index Kier, P.H., 255, 256, 268–271 Kinematics of elastic collision, 241 Knuth, D.E., 123, 124 Koebke, K., 190, 419 Kokotovic, P.V., 315 Koksma–Hlawka inequality, 137–141 Koksma, J.F., 137–141 Kolmogorov test, 125 Korobov, N.M., 130 Kramer, M.A., 316 Krylov methods, 2, 18, 63–75 biconjugate-gradient (BCG), 66 condition number, 67–69 conjugate-gradient, 67, 71, 72 generalized minimum-residual (GMRES), 66, 67, 70, 72 Krylov vectors, 64–66 minimum polynomial, 64 minimum-residual (MINRES), 66 preconditioner, 63, 68–71 preconditioning, 67–70 quasi-minimum residual (QMR), 66 Ritz–Galerkin approximation, 65, 66

L Lagrange multiplier, 107 Lamb, W.E., 232 Lane, M., 221 Laplace, M.P.-S., 117 Larsen, E.W., 54 Lathrop, K.D., 1 Latin Hypercube sampling, 302, 304, 305 Lattice structure, 126, 127 Lattice test, 127 Lawrence, R.D., 199, 204 Leal, L.C., 240 Lecot, C., 131 L’Ecuyer, P., 127, 128 Ledinegg, M., 443, 444 Lee, C.E., 171 Lee-Whiting, G.E., 237 Legendre addition theorem, 97 Legendre polynomials, 97, 102 Lehmer, D.H., 123 Lehner, J., 406, 408 Lehn, J., 127 Levine, M.M., 266 Levitt, L.B., 282, 283 Linear functional, 334, 341 Linear operators, 315, 317–319, 321 Local adjoint sensitivity analysis procedure, 291, 316, 317, 321–346, 349, 350 Local analysis, 302, 318, 320, 349

463 Lorentz, R., 419 Low discrepancy sequences, 120, 129–141, 149, 150, 152, 158. See also Quasi-random sequence good lattice points, 130, 140 Halton sequence, 130, 131, 136, 137 randomized, 141 Van der Corput sequence, 125, 130, 135–137 Lukyanov, A.A., 285 Lynn, J.E., 221

M Maillard, J.M., 283, 284 Maize, E.H., 149 Marchuk, G.I., 218 Markov chain Monte Carlo, 119 Markov property, 135, 137 Marleau, G., 215 Marsaglia, G., 123, 124, 127 Material temperature, 15, 16 Matrix problem, 144, 150, 155 symmetric, 96 Maynard, C.W., 147 MCNP, 22 Mean, 117, 123, 129, 134–136. See also Expected value free path, 154 sample, 120, 124, 133, 134, 136, 143, 155 value, 298, 339 Measurement equation, 296, 297 Measurements, 291–301, 349 direct, 297 equation, 296, 297 errors, 294–296 inaccuracy, 294 indirect, 295–297, 349 instruments, 293–295 theory, 293 Meromorphic, 225, 229, 230, 236 Methods brute-force, 314, 348 correlation ratio-based, 302, 308 fine mesh, 107 first-order reliability, 302, 307 sampling-based, 302–308 screening design, 302, 310–314, 346, 347 second-order reliability, 302, 307 Sobol’s, 302, 308, 309, 348 variance-based, 308–310 Mixed methods, 86, 87, 89

464 Model, 138 analytic, 155, 157 condensed history, 151, 153 deterministic, 155, 157 mathematical, 292, 295, 297, 307, 348 multiple scattering, 154 probability, 134, 155 problem, 117, 119, 137 Moldauer, P.A., 272 Moments, 118, 120, 153, 154, 157 of distribution, 297 equations, 157 higher-order, 297, 300 propagation equations, 298 Moment-to-discrete matrix, 41, 42 Monte Carlo, 20–23, 358, 360–361, 363–369, 371 Moore, M.S., 228–229, 237 Morokoff, W.J., 131 Morris, M.D., 312, 313 Moskowitz, B., 131 Multigroup, 1, 5–7, 12–14, 20, 49, 50, 61 cross section, 252–254 method, 85, 87, 96 Multiplication factor, 167, 170, 173 four factor formula, 181, 208 thermal utilization, 168 six factor formula, 168 Multistage procedure, 145

N Nayler, J., 265, 270 Nelkin, M., 232, 410 Nelkin model, 411 Neumann series, 118, 132 Neutral particle transport, 85 Neutron transport problem, 119 Newton’s method, 15 Niederreiter, H., 124 Nikolaev, M.N., 283 Nodal expansion method (NEM), 198 analytic nodal method, 199–200, 206, 207 nodal balance equation, 199, 204 Nodal method, 86, 106–110, 112 subelement, 112 variational, 106–110, 112 Nominal parameter value, 297, 314, 332, 340 Nonlinear functional, 301, 317, 324 Nonlinear operator, 317–319 Nordheim, L.W., 249, 250, 263–265, 267 Normal distribution, 134, 296 Normed vector space, 319 Nuclear chain reaction, 121

Index Nuclear reactor kinetics epilogue (1934–1999) adaptive-model reactor kinetics, 448–449 advanced reactors, 449–450 thermal hydraulic feedback, 447 twenty-first century, 450 exponential and non-exponential decay Bromwich contour, 415 effective multiplication constant, 411, 412 energy-spectrum regeneration mechanism, 412 Fredholm integral, 413 Green’s function, 413–414 Laplace transformed fission neutron production rate, 414 liquid metal fast breeder reactors (LMFBRs), 416, 417 pseudo mode and quasi-exponential die-away, 416 pseudo-mode decay, 412 Riemann surface, 415 slowing-down equation, 413–414 SUAK assemblies, 412, 413 time-dependent one-speed diffusion equation, 411 neutron thermalization, exponential decay, and diffusion cooling Boltzmann equation, 402 decay constant vs.buckling, 401 diffusion cooling phenomenon, 402 Maxwellian thermal neutron spectrum, 400, 401 neutron generator, 400 pulsed neutron experiment, 400 thermal reactor, 399 time asymptotic flux, 401 non-exponential decay and pulsed neutron die-away theory cross-section and time-of-flight measurements, 410 eigenvalue spectrum, 405 energy-dependent diffusion operator, 410 Fourier transformed Boltzmann equation, 407 fundamental space-angle mode, 406 image reactor theory, 407 Laplace transform, 405, 406 nonhomogeneous linear operator equation, 404 one-speed diffusion operator, 402, 403 point and continuous spectrum, 403

Index pseudo-exponential decay, 409 Riemann surface, 409 smack dead center, 411 spectrum components, 405 time-and velocity-dependent transport equation, 403 velocity-dependent Boltzmann equation, 408 point reactor kinetics equations credit and instantaneous neutrons, 380–381 Fermi and Wigner’s seminal contributions, 382 graphite/natural-uranium pile, 381 neutron and reactor physics, 378 number of radioactive nuclei, 379 one generation time, 380 one-group diffusion equation with delayed neutrons, 382–387 shape, time, and neutron importance, 388–393 variational and asymptotic formulations, 394–399 reactor dynamics boiling channels, 443 boiling water thermal reactors (BWRs), 443, 444 decay ratios (DRs), 445 density-wave oscillation, 444 Doppler effect, 441, 442 homogeneous equilibrium model, 444, 445 Hopf bifurcation, 444 Ledinegg instability, 443 nonlinear feedback mechanism, 441, 447 nonlinear reactor dynamics, 442, 443 poisoning effect, 442 power-flow map, 445, 446 space–time reactor kinetics adaptive-model kinetics calculations, 440–441 coarse-mesh and nodal methods, 429–434 equivalence theory, 419 finite-difference scheme, 419–423 homogenization theory, 434–440 steady-state diffusion equation, 417 synthesis methods, 418 time-dependent multigroup transport equation, 419 variational, modal, synthesis, and related methods, 423–429

465 time-dependent neutron diffusion equation elementary reactor physics, 376 hydrogenous materials, 377 law of diffusion, 377 neutron density, 378 via Panisperna Group, 378 Numerical Doppler-broadening heat equation approach, 240 kernel-broadening explicit, 238–239 implicit, 239–240 O Oil-well logging, 2, 14, 23, 29 borehole, 21, 22 gamma-gamma tool, 21, 22 neutron-neutron tool, 22 Okten, G., 131, 132, 140 Olson, A.P., 255, 271 One-at-a-time experiment, 302, 311, 312 Operator adjoint, 316, 330 linear, 315, 317–319, 321 nonlinear, 317–319 Organizing principles, 119–121 Osborn, R.K., 232 O’Shea, D.M., 235, 236 Outer iteration, 12, 17, 60–62 Owen, A., 131

P Parallel computing, 112 Parameters correlated, 301, 304, 306, 347, 348 nominal value, 297, 314, 332, 340 uncorrelated, 299–301, 304, 332 Pareto’s law, 311, 347 Parks, D.E., 232 Partial current, 110 Particle transport, 117, 119, 133 Peierls, R., 221, 223–225, 227–229 Perez, R.B., 228–230 Periodic function, 130, 140 Perturbation Monte Carlo, 151, 158 Perturbation theory, 315–317 Phase space, 130, 138, 142, 145, 146, 148–150, 156 Pisot, C., 125 Placzek, G., 243–246 Planck function, 15, 16 Planck’s constant, 15 Pn or PN method, 96, 102–104

466 Point reactor kinetics equations bare homogeneous reactor bi-orthogonality properties, 387 buckling mode, 383, 384, 386, 387 delayed-neutron precursor concentration, 382–383 fundamental spatial mode, 386 Helmholtz equation, 384 in-hour equation, 386 nonhomogeneous equation, 387 ordinary differential equations (ODEs), 385 partial differential equations (PDEs), 383 separation constant, 384 space-and time-dependent neutron number density, 382 credit and instantaneous neutrons, 380–381 Fermi and Wigner’s seminal contributions, 382 graphite/natural-uranium pile, 381 neutron and reactor physics, 378 number of radioactive nuclei, 379 one generation time, 380 shape, time, and neutron importance adiabatic method, 393 analogous decomposition, 389 Boltzmann equation, 388, 390 effective delayed neutron fraction, 392 instantaneous fundamental-mode eigenfunction, 393 lambda-and omega-mode, 392–393 normalization condition, 390–391 precursor equation, 390 reactivity, 391 shape and time function, 388 variational and asymptotic formulations asymptotic expansion technique, 395 Boltzmann equation, 394–395 energy-/speed-dependent transport theory, 399 Euler–Lagrange equation, 395 fundamental lambda-mode eigenfunction, 397 fundamental omega-mode eigenfunction, 398 generalized point reactor equations, 398 homogeneous transport equation, 396 instantaneous critical reactor, 397 neutron number density, 396 normalization condition, 394, 397 Polynomials complete, 108

Index piecewise, 105 trial functions, 94, 105 Pomeranchouk, I.Y., 218, 262 Pontecorvo, B., 377 Porter, C.E., 274 Practical representations Adler–Adler approximation, 226–228 generalized pole representation, 229–230 multilevel Breit–Wigner approximation, 225–226 Reich–Moore approximation, 228–229 single level Breit–Wigner approximation, 224–225 Primal method, 89–91 Probability density, 144–147, 149, 151, 153, 154, 156 conditional, 135 Propagation of errors, 291, 293, 297–301, 318, 332, 349 Pseudorandom number, 120, 124, 126. See also Random number Purely random errors, 295, 296

Q Quadrature, numerical, 101, 111 Quantile, 304–306 Quasi-Monte Carlo method, 129, 130, 134–140, 142, 146, 149, 155 Quasi-random, 124, 131, 134–142, 157 errors, 296 sequence, 120, 125, 129, 133–135, 137, 139 (see also Low discrepancy sequence)

R Radiation, thermal, 112 Radical inverse function, 130 Random errors, 294–296 Random number, 120, 122, 123, 128, 152. See also Pseudorandom number periodic sequences, 125 Random sampling, 302, 304, 305 Random variable, 118, 123, 125, 132, 134, 135, 142, 143, 152, 156. See also Estimator Random walk, 117–119, 130, 132, 133, 135–138, 140, 144, 146–152, 155 Rasetti, F., 377 Ray effects, 2, 104 Rayleigh–Ritz method, 424, 427 Ray tracing, 111 Reactor core database, 204, 205

Index Rebalance, coarse-mesh rebalance, 13 Reciprocity, 147 Reduced source, 145, 146 Reich, C.W., 228–229, 237 Relative error, 134, 136 Reproducing kernel Hilbert space, 137 Resonance absorption in heterogeneous media closely-packed lattice cylindrical configuration, 265 Dancoff effect, 263, 265 infinite slabs, 264–265 Nordheim’s approximation, 264 rational approximation, 265–266 two region cell Chernick’s method, 257 general feature of collision probability, 258–260 rational approximation-based method NR-approximation, 260–262 NRIM-approximation, 262 Wigner’s rational approximation, 259–261, 264 Resonance escape probability, 143 Resonance integral approximations in homogeneous media infinite mass approximation (NRIM), 248 intermediate resonance (IR) approximation, 248–249 narrow resonance (NR) approximation, 247–248 Nordheim’s method, 249–250 various J-integrals J(™k ; “k ; ak ; bk ), 251, 252, 280 J(™r ; “r ), 248, 262 J(™r ; “r ; ar ), 247, 280 Response, 292, 293, 296–298, 300–303, 305–308, 310–321, 323–329, 331–333, 336, 337, 339, 340, 342–344, 346–349 Ribon, P., 283, 284 Ricatti equation, 339–342 Rigorous resonance absorption treatment cylindrical unit cells, 268–270 homogeneous media, 255 unit cell with many plates, 271–272 unit cell with multiple regions, general, 268 Ritz–Galerkin method, 425 R-matrix theory Kapur–Peierls version, 224–225 Wigner–Eisenbud version channel matrix, 221 hard-sphere phase shift factor, 222 level matrix, 221, 222 level shift factor, 222, 223

467 penetration factor, 222 spherical Hankel function, 222 Robba, A.A., 255, 256, 268–271 Ronen, Y., 349 Rosenfeld, L., 230 Rosenwasser, E., 315 Rothenstein, W., 263–265, 271 Roussopoulos, P., 425, 426, 449 Roy, R., 215 Russian roulette, 147, 148

S Saltelli, A., 309 Sample space, 146, 149 Sampling Latin Hypercube, 302, 304, 305 random, 302, 304, 305 stratified importance, 302 Sampling-based methods, 302–308 Sandwich rule, 299, 300, 332, 337 Saunderson, J.L., 154 Scalar flux, 87, 88, 107, 110, 111 Scattering anisotropic, 86, 87, 96–99 isotropic, 86, 87, 110 matrix, 40, 41 Scattering operator, 62 anisotropic scattering, 2, 20, 30, 36, 38–45, 59, 70 Boltzmann scattering, 19, 39, 43, 44, 48, 49 continuous scattering, 45–49 continuous-slowing-down, 19, 20, 45, 49–51 Thompson scattering, 15 upscattering, 12 Screening design method, 302, 310–314, 346, 347 Second-order reliability method, 302, 307 Segrè, E., 377 Self-shielding effects degree of, 220 self-shielding factor, 220 Sensitivity analysis deterministic, 291, 302, 314–330, 347 statistical, 291, 302–314, 318, 346, 347 Sensitivity coefficient, 306, 347 Sequential bifurcation design, 302, 313 Sequential Monte Carlo, 144. See also Adaptive Monte Carlo correlated sampling, 145, 148, 151 Serial correlation, 123, 125, 127 Simplified angular approximation, 102–104, 107

468 Sirakov, I.A., 285 Slowing-down in homogeneous media Placzek oscillations, 244–246 slowing-down density, 243–244 slowing-down equation, 258, 260, 267 Smith, K.S., 419 Sn or SN , 96, 102–104 Sobol’s method, 302, 308, 309, 348 Solbrig, A.W., 231, 233–235, 238, 239 Source external, 96 group, 98 Source iteration, 13, 19, 45, 49, 51, 53–55, 57, 60, 61, 63 marching, 11 transport sweep, 11, 12 update equation, 11, 58 Space–time reactor kinetics adaptive-model kinetics calculations, 440–441 coarse-mesh and nodal methods accurate homogenization method, 433 constraint equations, 432 discrete-variable unknowns, 429 higher-order nodal method, 433, 434 hybrid nodal Green’s function method, 432 neutron diffusion and transport equation, 434 ODE, 431 one-group diffusion equation, 430 partial differential equation, 431 quadratic transverse surface fit, 431, 432 equivalence theory, 419 finite-difference scheme coupled delayed neutron precursor equation, 420 forward and backward Euler scheme, 419 g-th group diffusion removal equation, 421 light water reactor (LWR), 422, 423 linear algebraic equation, 421 multigroup diffusion equations, 420 theta-difference method, 422 homogenization theory direction-dependent diffusion coefficients, 434 exact equivalence theory, 434, 435 flux discontinuity factors, 435, 436 forward and adjoint shape functions, 439 fuel-assembly diffusion equation, 437

Index generalized equivalence theory, 435 lattice-cell transport equation, 436 nodal Green’s method, 438 self-consistent models, 439 variable-time-step adaptive-grid calculation, 440 steady-state diffusion equation, 417 synthesis methods, 418 time-dependent multigroup transport equation, 419 variational, modal, synthesis methods Euler–Lagrange equation, 424 forward and adjoint equation, 426 linear function space approximation method, 427 linear nonhomogeneous algebraic system, 426 linear self-adjoint equation, 423, 424 multichannel synthesis methods, 428, 429 multigroup diffusion and transport equations, 425 one-group steady-state diffusion equation, 423 Roussopoulos function, 425, 426 trial function, 424 variational synthesis methods, 428 weight functions, 427 Spanier, J., 131, 132 Spherical harmonics, 86, 96, 100, 102–104, 107, 154 odd parity, 97 simplified, 103 Splitting, 147, 148 SPn, 96, 102, 104, 108 SSn, 96, 102, 103, 108 Stacey, W.M., Jr., 254 Standard deviation, 120, 134, 136, 142 sample, 134, 140 State of the art, 140–141 for error reduction, 154 for foundations/theoretical developments, 121, 157 for generating sequences, 133 Stationary condition, 91, 95, 100, 101, 106, 107, 109, 111 Statistical distributions level width Dyson’s distribution, 274, 275 Dyson’s two-level correlation function, 276 Gauss orthogonal ensemble, 275

Index level correlation function, 275–276 Wigner distribution, 274 partial width, 274 Statistical error, 141, 155 Statistical independence, 137 Statistical test, 120, 123, 125, 134 Statistical theory fundamentals addition theorem, 273 basic rule, 273 multiplication theorem, 273 Steen, N.M., 237 Stratified importance sampling, 302 Stratified sampling, 142, 149–150 Subcritical, 359, 365, 366, 368 limit, 355 Subjective uncertainty, 305 Superposition principle, 143, 144 Survival biasing, 150 System parameter, 298–300, 320, 330, 332, 348 response, 317, 319 Systematic errors, 294, 296 Systematic fractional replicate design, 302, 312 Systematic sampling, 118 Systematic source sampling, 149

T Takahashi, H., 265 Tangent linear model, 320, 330 Taylor series, 297, 298, 350 Test function, 88 Thacher, H.C., 235, 236 Theoretical foundations, 119–121, 155–157 Thermal radiation transport, 14–18, 35 Thomas, R.G., 221, 274 Tomovic, R., 315 Transmission probability, 259, 263, 265 Transport, electron, 112 Transport equation, 118, 135, 141, 142, 144–149, 151, 153, 155–157 adjoint, 146, 148, 156, 157 Transverse-integration nodal method, 429–431 Treatment of unresolved resonances average unshielded cross sections, 278–279 conceptual aspects, 276–277 fluctuation integral, 278–279 maximization of information entropy, 286 probability table method characteristic function approach, 286 conceptual basis, 282–283

469 integral transform technique, 286 moment method, 284–285 resonance ladder method, 283–284 self-shielded cross sections, 276, 279–281 strength function, 278 Trial function, 94, 95, 110, 112 spatial, 105 Truncation error, 5, 13, 30, 59 Turing, A.M., 235, 236

U Ultra-fine group, 253 Uncertainty, 291–350 epistemic, 293 measurement, 293–297 stochastic, 292, 293 subjective, 305 Uncorrelated parameters, 299–301, 304, 332 Uniform sequence, 122, 140 Unit cell, 174, 175, 191, 192, 206 Unresolved boundary layers, 14, 20, 35 Use of expected values, 142, 150–151 Ussachoff, L.N., 388, 393, 394, 397

V Vacuum regions, 112 Van der Corput, J.G., 125, 135–137 Variance, 118, 140, 143, 144, 147, 148, 155, 157, 298, 299, 302, 305, 306, 308–310, 312, 331–339 of population, 134 of the variance, 134 Variance-based method, 308–310 Variance reduction, 141–143, 147–150 adaptive importance sampling, 148, 158 antithetic variates, 143–144 control variates, 143–146 exponential transformation, 147 importance sampling, 138, 141, 142, 145–149, 154, 155, 158 sequential correlated sampling, 145, 148, 158 Variational method, 85–87, 91–93, 99, 103, 106, 107 Various J-integrals evaluations, Gauss–Chebyshev quadrature, 251 Jk *, 251 Vladimirov, V.S., 85, 91 Vucobratovic, M., 315

470 W Wang, J., 131 Warsa, J.S., 75 Weak form, 86–93, 99 Weighted residuals, 90, 106 Weight windows, 148, 151 Weinberg, A.M., 379, 380 Whitesides, G.E., 364, 371 Wigner, E.P., 218, 219, 221–225, 230, 232, 246–248, 250, 251, 259, 261–263,

Index 265, 267, 272, 274–276, 316, 318, 375, 376, 378–382 Williams, M.M.R., 168 Wing, G.M., 406, 408 Winkler, R.L., 304

Y Yip, S., 232 Yusupov, R., 315