Duderstadt Hamilton Nuclear Reactor Analysis
Duderstadt Hamilton Nuclear Reactor Analysis
Reactor
Analysis
James J. Duderstadt
Louis J. Hamilton
Department of Nuclear Engineering
The Uniwrsity of Michigan
Ann A r h r , Michigan
James J. Duderstad1
Contents
Louis J. Hamilton
PART 1
Introductory Concepts
of Nuclear Reactor
Analysis
I. Nuclear Reactions
11. Nuclear Fission
xiii
xiv
An Introduction
to Nuclear Reactor
PART 3 Core Design
The Multigroup
Diffusion Method CHAPTER 11 General Aspects of Nuclear Reactor Core Design 447
I. Introduction
CHAPTER 14 Reactivity Control 11. Properties of the Dirac &Function
A. Alternative Representations
I. Introduction B. Properties
11. Movable Control Rods C. Derivatives
111. Burnable Poisons
IV. Chemical Shim D. Some Properties of Special Functions
V. Inherent Reactivity Effects
E. Some Assorted Facts on Linear Operators
CHAPTER 15 Analysis of Core Composition Changes
I. Scalar Products
I. Fission Product Poisoning 111. Linear Operators
11. Fuel Depletion Analysis 111. Linear Vector Spaces
111. Nuclear Fuel Management IV. Properties of Operators
V. Differential Operators
I. Some Definitions
11. Matrix Algebra
I. Motivation
11. "Cookbook" Laplace Transforms
It has been more than three decades since the first nuclear reactor achieved a
critical fission chain reaction beneath the old Stagg Field football stadium at the
University of Chicago. Since that time an extensive worldwide effort has been
directed toward nuclear reactor research and development in an attempt to harness
the enormous energy contained within the atomic nucleus for the peaceful genera-
tion of power. Nuclear reactors have evolved from an embryonic research tool into
the mammoth electrical generating units that drive hundreds of central-station
power plants around the world today. The recent shortage of fossil fuels has made
it quite apparent that nuclear fission reactors will play a dominant role in meeting
man's energy requirements for decades to come.
For some time electrical utilities have been ordering and installing nuclear plants
in preference to fossil-fueled units. Such plants are truly enormous in size, typically
generating over 1000 MWe (megawatts-electric) of electrical power (enough to
supply the electrical power needs of a city of 400,000 people) and costing more
than one billion dollars. It is anticipated that some 500 nuclear power plants
will be installed in the United States alone by the year 2000 with an electrical
generating capacity of about 500,000 MWe and a capital investment of more than
$600 billion,' with this pattern being repeated throughout the world. The motiva-
tion for such a staggering commitment to nuclear power involves a number of
factors that include not only the very significant economic and operational advan-
tages exhibited by nuclear plants over conventional sources of power, but their --
substantially lower environmental impact and vastly larger fuel resources as well.L-'
The dominant role played by nuclear fission reactors in the generation of
electrical power can be expected to continue well into the next century. Until at
least A.D. 2000, nuclear power will represent the only viable alternative to
4 / INTRODUCTORY CONCEPlS OF NUCLEAR POWER REACTOR ANALYSIS INTRODUCTION TO NUCLWR POWER GENERATION / 5
fossil-fueled plants for most nation^.^.^ The very rapid increase in fossil-fuel costs The radiation produced by reactors can also be used to transmute nuclei into
that has accompanied their dwindling reserves has led to a pronounced cost artificial isotopes that can then be used, for example, as radioactive tracers in
advantage for nuclear plants which is expected to widen even further during the industrial or medical applications. Reactors can use the same scheme to produce
next few Of course, there are other longer-range alternatives involving nuclear fuel from nonfissile materials. For example, 238U can be irradiated by
advanced technology such as solar power, geothermal power, and controlled neutrons in a reactor and transmuted into the nuclear fuel 2 3 9 ~ This
~ . is the process
thermonuclear fusion. However the massive practical implementation of such utilized to "breed" fuel in the fast breeder reactors currently being developed for
alternatives, if proven feasible, could probably not occur until after the turn of the commercial application in the next decade.
century since experience has shown that it takes several decades to shift the energy Small, compact reactors have been used for propulsion in submarines, ships,
industry from one type of fuel to another? due both to the long operating lifetime aircraft, and rocket vehicles. Indeed the present generation of light water reactors
of existing power machinery and the long lead times needed to redirect manufac- used in nuclear power plants are little more than the very big younger brothers of
turing capability. Hence nuclear fission power will probably be the dominant new the propulsion reactors used in nuclear submarines. Reactors can also be utilized as
source of electrical power during the productive lifetimes of the present generation small, compact sources of long-term power, such as in remote polar research
of engineering students. stations or in orbiting satellites.
Yet by far the most significant application of nuclear fission reactors is in large,
central station power plants. A nuclear power plant is actually very similar to a
I. NUCLEAR FISSION REACTORS fossil-fueled power plant, except that it replaces the coal or oil-fired boiler by a
nuclear reactor, which generates heat by sustaining a fission chain reaction in a
The term nuclear reactor will be used in this text to refer to devices in which suitable lattice of fuel material. Of course, there are some dramatic differences
controlled nuclear fission chain reactions can be maintained. (This restricted between a nuclear reactor and, say, a coal-fired boiler. However the useful quantity
definition may offend that segment of the nuclear community involved in nuclear produced by each is high temperature, high pressure steam that can then be used to
fusion research, but since even a prototype nuclear fusion reactor seems several run turbogenerators and produce electricity. At the center of a modern nuclear
years down the road, no confusion should result.) In such a device, neutrons are plant is the nuclear steam supply system (NSSS), composed of the nuclear reactor,
used to induce nuclear fission reactions in heavy nuclei. These nuclei fission into its associated coolant piping and pumps, and the heat exchangers ("steam genera-
lighter nuclei (fission products), accompanied by the release of energy (some 200 tors") in which water is turned into steam. (A further glance at the illustrations in
MeV per event) plus several additional neutrons. These fission neutrons can then be Chapter 3 will provide the reader with some idea of these components.) The
u t i l i to induce still further fission reactions, thereby inducing a chain of fission remainder of the power plant is rather conventional.
events. In a very narrow sense then, a nuclear reactor is simply a sufficiently large Yet we must not let the apparent similarities between nuclear and fossil-fueled
mass of appropriately fissile material (e.g., 235Uor u9Pu) in which such a controlled power plants overshadow the very significant differences between the two systems.
fission chain reaction can be sustained. Indeed a small sphere of "5U metal slightly For example, in a nuclear plant sufficient fuel must be inserted into the reactor
over 8 cm in radius could support such a chain reaction and hence would be core to allow operation for very long periods of time (typically one year). The
classified as a nuclear reactor. nuclear fuel cycle itself is extremely complex, involving fuel refining, fabrication,
However a modern power reactor is a considerably more complex beast. It must reprocessing after utilization in the reactor, and eventually the disposal of radioac-
not only contain a lattice of very carefully refined and fabricated nuclear fuel, but tive fuel wastes. The safety aspects of nuclear plants are also quite different, since
must as well provide for cooling this fuel during the course of the chain reaction as one must be concerned with avoiding possible radiological hazards. Furthermore
fission energy is released, while maintaining the fuel in a very precise geometrical the licensing required by a nuclear plant before construction or operation demands
arrangement with appropriate structural materials. Furthermore some mechanism a level of sophisticated analysis totally alien to fossil-fueled plant design.
must be provided to control the chain reaction, shield the surroundings of the Therefore even though the NSSS contributes only a relatively modest fraction of
reactor from the intense nuclear radiation generated during the fission reactions, the total capital cost of a nuclear power plant (presently about 20%), it is of central
and provide for replacing nuclear fuel assemblies when the fission chain reaction concern since it not only dictates the detailed design of the remainder of the plant,
has depleted their concentration of fissile nuclei. If the reactor is to produce power but also the procedures required in plant construction and operation. Furthermore
in a useful fashion, it must also be designed to operate both ewnomically and it is the low fuel costs of the NSSS that are responsible for the economic
safely. Such engineering constraints render the actual nuclear configuration quite advantages presently enjoyed by nuclear power generation.
complex indeed (as a quick glance ahead to the illustrations in Chapter 3 will The principal component of the NSSS is, of wurse, the nuclear reactor itself. A
indicate). rather wide variety of nuclear reactors are in operation today or have been
Nuclear reactors have been used for over 30 years in a variety of applications. proposed for future development. Reactor types can be characterized by a number
They are particularly valuable tools for nuclear research since they produce of features. One usually distinguishes between those reactors whose chain reactions
copious amounts of nuclear radiation, primarily in the form of neutrons and are maintained by neutrons with characteristic energies comparable to the energy
gamma rays. Such radiation can be used to probe the microscopic structure and of thermal vibration of the atoms comprising the reactor core (thermal reactors)
dynamics of matter (neutron or gamma spectroscopy). and reactors in which the average neutron energy is more characteristic of the
6 / JNTRODU(T0RY CONCEPTS OF NUCLEAR POWER REACTOR ANALYSIS INTRODUCllON TO NUCLEAR POWER GENERATION / 7
much higher energy neutrons released in a nuclear fission reaction (fast reactors). assist in both the nuclear design of fission reactors and their integration into large
Yet another common distinction refers to the type of coolant used in the reactor. power systems. In the early days of the reactor industry a nuclear engineer was
In the United States, and indeed throughout the world, the most popular of the usually regarded as a Ph. D.-level reactor physicist primarily concerned with
present generation of reactors, the light water reactor (LWR) uses ordinary water nuclear reactor core research and design. Today, however, nuclear engineers are
as a coolant. Such reactors operate at very high pressures (approximately 70-150 needed not only by research laboratories and reactor manufacturers to develop and
bar) in order to achieve high operating temperatures while maintaining the water in design nuclear reactors, but also by the electrical utilities who buy and operate the
its liquid phase. If the water is allowed to boil in the core, the reactor is referred to nuclear power plants, and by the engineering companies who build the power
plants and service them during their operating lifetimes.
as a boiling water reactor (BWR), while if the system pressure is kept sufficiently
Hence an understanding of core physics is not sufficient for today's nuclear
high to prevent bulk boiling (155 bar), the reactor is known as a pressurized water
reactor (PWR). Such reactors have benefited from a well-developed technology and engineer. He must also learn how to interface his specialized knowledge of nuclear
performance experience achieved in the nuclear submarine program. reactor theory with the myriad of other engineering demands made upon a nuclear
A very similar type of reactor uses heavy water ( 4 0 ) either under high pressure power reactor and with a variety of other disciplines, including mechanical,
as a primary coolant or simply to facilitate the fission chain reaction. This electrical, and civil engineering, metallurgy, and even economics (and politics), just
particular concept has certain nuclear advantages that allow it to utilize low- as specialists of these other disciplines must learn to interact with nuclear en-
enrichment uranium fuels (including natural uranium). It is being developed in gineers. In this sense, he must recognize that the nuclear analysis of a reactor is
Canada in the CANDU series of power reactors and in the United Kingdom as only one facet to be considered in nuclear power engineering. To study and master
steam generating heavy water reactors (SGHWR). it outside of the context of these other disciplines would be highly inadvisable. In
Power reactors can also utilize gases as coolants. For example, the early MAG- the same sense, those electrical, mechanical, or structural engineers who find
NOX reactors developed in the United Kingdom used low-pressure CO, as a themselves involved in various aspects of nuclear power station design (as ever
coolant. A particularly attractive recent design is the high-temperature gas-cooled increasing numbers are) will also find some knowledge of nuclear reactor theory
reactor (HTGR) manufactured in the United States which uses high-pressure useful in the understanding of nuclear components and interfacing with nuclear
helium. Related gas-cooled reactors include the pebble-bed concept and the ad- design.
vanced gas cooled reactors (AGR) under development in Germany and the United Future nuclear engineers must face and solve complex problems such as those
Kingdom, respectively. involved in nuclear reactor safety, environmental impact assessment, nuclear power
All of the above reactor types can be classified as thermal reactors since their plant reliability, and the nuclear fuel cycle, which span an enormous range of
fission chain reactions are maintained by low-energy neutrons. Such reactors disciplines. They must always be concerned with the economic design, construc-
comprise most of the world's nuclear generating capacity today, and of these the tion, and operation of nuclear plants consistent with safety and environmental
LWR is most common. It is generally agreed that the LWR will continue to constraints. An increasing number of nuclear engineers will find themselves con-
dominate the nuclear power industry until well into the 1980s, although its market cerned with activities such as quality assurance and component standardization as
may tend to be eroded somewhat by the successful development of the HTGR the nuclear industry continues to grow and mature, and of course all of these
or advanced heavy water reactors. problems must be confronted and handled in the public arena.
However as we will see in the next chapter, there is strong incentive to develop a
fast reactor which will breed new fuel while producing power, thereby greatly 111. SCOPE OF THE TEXT
reducing nuclear fuel costs. Such fast breeder reactors may be cooled by either
liquid metals [the liquid metal-cooled fast breeder reactor (LMFBR)] or by helium Our goal in this text is to develop in detail the underlying theory of nuclear
[the gas-cooled fast breeder reactor (GCFR)]. Although fast breeder reactors are fission reactors in a manner accessible to both prospective nuclear engineering
not expected to make an appreciable impact on the nuclear power generation students and those engineers from other disciplines who wish to gain some
market until after 1990, their development is actively being pursued throughout the exposure to nuclear reactor engineering. In every instance we attempt to begin with
world today. the fundamental scientific principles governing nuclear fission chain reactions and
Numerous other types of reactors have been proposed and studied-some even then carry these fundamental concepts through to the level of realistic engineering
involving such exotic concepts as liquid or gaseous fuels. Although much of the applications in nuclear reactor design. During this development we continually
analysis presented in this text is applicable to such reactors, our dominant concern stress the interplay between the nuclear analysis of a reactor core and the parallel
is with the solid-fuel reactors cooled by either water, sodium, or helium, since these nonnuclear design considerations that must accompany it in any realistic nuclear
will comprise the vast majority of the power reactors installed during the next reactor analysis.
several decades. We must admit a certain preoccupation with nuclear power reactors simply
because most nuclear engineers will find themselves involved in the nuclear power
11. ROLE OF THE NUCLEAR ENClNEER industry. This will be particularly apparent in the examples we have chosen to
discuss and the problems we have emphasized. However since our concern is
The nuclear engineer will play a very central role in the development and always with emphasis of fundamental concepts over specific applications, most of
application of nuclear energy since he is uniquely characterized by his ability to the topics we develop have a much broader range of validity and would apply
8 / INTRODUCTORY CONCEPTS OF NUCLEAR POWER REACTOR ANALYSIS
INTRODUCTION TO NUCLEAR POWER GENERATION / 9
equally well to the analysis of other types of nuclear reactors. And although our 4. A study of base-load alternatives for the Northeast Utilities System, Report to the Board
principal target is the prospective nuclear engineer, we would hope that engineers of Trustees of Northeast Utilities, Arthur D. Little, Inc. (1973).
from other disciplines would also find this text useful as an introduction to the 5. Nuclear Fuel Resources and Requirements, USAEC Report WASH-1243 (1973); Nucl.
concepts involved in nuclear reactor analysis. Ind. 21 (2) (1974); Nucl. News 17, (5) (1974).
The present text develops in four progressive stages. Part 1 presents a very brief 6. David J. Rose, Science 184, 351 (1974).
introduction to those concepts from nuclear physics relevant to nuclear fission 7. An Assessment of Accident Risks in U. S. Commercial Nuclear Power Plants, USAEC
reactors. These topics include not only a consideration of the nuclear fission Report WASH-1400 (1974).
process itself, but also a consideration of the various ways in which neutrons, which Most current information concerning the nuclear power industry appears in a number of
act as the carrier of the chain reaction, interact with nuclei in the reactor core. We journal publications. Although the number of journals appearing in the field of nuclear
next consider from a qualitative viewpoint the general concepts involved in science and engineering is quite voluminous, a brief list of journals of more general interest
studying nuclear chain reactions. Part 1 ends with an overview of nuclear reactor would include:
engineering, including a consideration of the various types of modern nuclear
Nuclear Engineering International (Europressatom): A British journal distinguished by its
reactors, their principal components, and a qualitative discussion of nuclear reactor elaborate, multicolored diagrams of nuclear power plants.
design. Nuclear Industry (Atomic Industrial Forum): A monthly news magazine written more from
Parts 2-4 are intended to develop the fundamental scientific principles underly- the viewpoint of the consumer of nuclear power products-namely, the electrical utilities.
ing nuclear reactor analysis and to apply these principles for derivation of the most Nuclear News (American Nuclear Society): A monthly news magazine published by the
common analytical tools used in contemporary reactor design. By way of illustra- American Nuclear Society, the principal technical organization of nuclear engineering in
tion, these tools are then applied to analyze several of the more common and the United States.
significant problems facing nuclear engineers. Nuclear Safety (Office of Information Services, U. S. Atomic Energy Commission): A
Part 2 develops the mathematical theory of neutron transport in a reactor. It journal highlighting recent developments in the field of nuclear reactor safety.
begins with the most general description based on the neutron transport equation Nuclear Science and Engineering (American Nuclear Society): The principal technical
and briefly (and wry qualitatively) reviews the standard approximations to this research journal of nuclear engineering.
equation. After this brief discussion, we turn quickly to the development of the
simplest nontrivial model of a nuclear fission reactor, that based upon one-speed Other more research-oriented journals in the field of nuclear science and engineering
include:
neutron diffusion theory. This model is used to analyze both the steady state and
time-dependent behavior of nuclear reactors, since although the model has very Annals of Nuclear Science and Engineering (formerly Journal of Nuclear Energy) (Pergamon,
limited validity in practical reactor analysis, it does illustrate most of the concepts New York).
as well as the calculational techniques used in actual reactor design. Journal of Nuclear Science and Technologv (Atomic Energy Society of Japan).
In Part 3, we develop the principal tool of modern nuclear reactor design, the Nuclear Engineering and h i g n (North-Holland, Amsterdam).
multigroup diffusion model. Particular attention is devoted to the calculation of the Nuclear Technology (American Nuclear Society).
Soviet Atomic Energy (Consultant's Bureau).
multigroup constants appearing in these equations, as well as to the practical
numerical solution of the equations themselves. References on nuclear reactor theory include:
In Part 4, we attempt to give an overview of the methods used in nuclear reactor
core design. In particular, we consider the application of the concepts and tools Bell, G. I., and Glasstone, S., Nuclear Reactor Theory, Reinhold, New York (1970).
developed in the earlier sections to a variety of problems faced by the nuclear Henry, Allan F., Nuclear Reactor Analysis, M.I.T. Press, Cambridge, Mass. (1975).
engineer, including criticality calculations, the determination of core power distri- Glasstone, S. and Edlund, M. C., The Elements of Nuclear Reacror Theory, Reinhold, New
York (1952).
butions and thermal-hydraulics analysis, burnup and control studies, and fuel-
Glasstone, S. and Sesonske, A., Nuclear Reactor Engineering, 2nd Edition. Van Nostrand,
loading requirements. While certainly incomplete, we do feel that the problems we Princeton, N.J.(1975).
have chosen to examine are representative of those encountered in nuclear reactor Lamarsh, J. R., Introduction to Nuclear Reactor Theoty, Addison-Wesley, Reading, Mass.
design and serve to illustrate the concepts developed in the earlier chapters of the (1966).
text. Lamarsh, John R., Introduction lo Nuclear Engineering, Addison-Wesley, Reading, Mass.
(1975).
Meghreblian. R. V. and Holmes, D. K., Reactor Analysis, McGraw-Hill, New York (1960).
REFERENCES Murray, Raymond L., Nuclear Energv, Pergamon Press, Nevr York (1975).
- -
Weinberg, A. M. and Wigner, E. P., The Physical Theory of Neutron Chain Reactors, Chicago
1. The Nuclear Industry, USAEC Report WASH-1 17473, 1973. The USAEC publishes an U. P. (1958).
annual report on the status of private nuclear industry within the United States. These Zweifel, P. F., Reactor Physics, McGraw-Hill, New York (1973).
reports provide a very comprehensive survey of the growth of the nuclear power industry.
2. Chauncey Starr, Sci. Amer. 225, 39 (1971).
3. Fourth International Conference on the Peaceful Uses of Atomic Energy (International
Atomic Energy Agency, Geneva, 1972). Vol. 1.
THE NUCLEAR PHYSICS OF FISSION CHAIN REACTIONS / I1
propagate a c h i n of fission reactions. In this sense, then, the neutron plays the role
The Nuclear Physics of of the chain carrier while the fission reactions supply the desired energy.
However there are other possible nuclear reactions a neutron can undergo that
Fission Chain Reactions do not lead to fission and hence are unproductive in nature. Indeed since there are
usually two or three neutrons emitted in each fission reaction, it should be apparent
that if each neutron resulted in another fission, the chain reaction would quickly
grow without bound. One such parasitic reaction involves the capture of the
neutron by a nucleus which then emits a gamma ray rather than fissioning.
Another possible reaction involves the neutron simply bouncing or scattering off of
a nucleus. After several such scattering reactions, the neutron might eventually leak
out of the uranium core of the reactor. Such processes remove neutrons from the
reactor and tend to inhibit the chain reaction.
Therefore one of the primary tasks of the nuclear engineer is to follow the
neutron "economy" in a nuclear reactor in order to monitor and control the
behavior of the fission chain reaction. That is, he must learn how to design the
reactor so that there is a balance between the production of neutrons in fission
reactions and the loss of neutrons due to capture or leakage. The study of such
processes is known as either nuclear reactor theory, nuclear reactor physics or
sometimes simply as neutronics. It is essentially the subject of this text.
However the achievement of a stable chain of fission reactions is only a part of
the responsibility of the nuclear engineer. In addition he must learn how to extract
and use the energy liberated in these fission reactions. This task involves the
subjects of heat transfer, fluid flow, structural and materials analysis, and power
systems analysis and interacts strongly with the nuclear analysis of a reactor a r e .
The primary objective in the design and operation of a nuclear reactor is the It is discussed in the latter chapters of the text.
utilization of the energy or radiation released by a controlled chain reaction of
We first turn our attention to a development of the fundamental concepts
nuclear fusion events maintained within the reactor core. Such fission reactions
involved in predicting the distribution of neutrons in a nuclear reactor in order to
occur when a heavy atomic nucleus such as ='U splits or fissions into two lighter
nuclei with an attendant release of both energy and radiation. Yet just how are understand and design a fission chain reaction system. We need to consider
such fission reactions induced in a reactor? The rate at which such naturally essentially two different subjects: (a) the determination of the probabilities of
occurring heavy nuclei will fission on their own (spontaneousfission) is very slow. It occurrence of various neutron-nuclear reactions and @) the derivation and solu-
should also be apparent that one cannot simply smash two nuclei together to tion of an equation that uses these probabilities to determine the neutron density
induce such a reaction, since the large electrical charges of heavy nuclei would lead and fission reaction rate in a nuclear reactor core.
The above discussion clearly indicates the importance of being able to determine
to a very strong repulsion. A more attractive idea is to slam a neutral particle
the rate at which various types of neutron-nuclear reactions occur within the
(which doesn't feel the nuclear charge) into a big, "overweight" nucleus and hope reactor. However it is important to keep in mind that there are enormous numbers
that this splits it. An ideal candidate for the incident particle is the neutron. Indeed of neutrons (typically l d per cm3) and even larger numbers of nuclei ( I d 2 per cm3)
experiments have shown that certain nuclei have an enormous appetite for in the reactor core. Hence we really need concern ourselves only with the average
neutrons, but after swallowing them suffer a case of violent indigestion that results behavior of the neutrons and nuclei in the reactor in a statistical sense. That is, we
in their fission. As an example of such a reaction, consider a neutron incident upon wish to calculate the probabilities that various types of neutron-nuclear interactions
a ='U nucleus: will occur. These reaction probabilities are expressed in terms of parameters called
+
Neurron 23S~+fission
products + more neutrons + energy. nuclear cross sections.
These cross sections represent the fundamental data utilized by the nuclear
The products of such a reaction (e.g., lighter nuclei, neutrons, and gammas) emerge engineer in his analysis of a nuclear reactor, much in the same way that thermal or
with very large kinetic energy (some 200 MeV) which is then converted into heat as structural data are used by the mechanical engineer or circuit device parameters
they slow down by banging into neighboring atoms in the reactor fuel. It is this are used by the electrical engineer. Hence some familiarity with the physics
heat energy that one utilizes to produce steam and eventually electrical power in a underlying the determination and behavior of such cross sections is necessary for
nuclear power plant. effective nuclear reactor analysis.
Yet just as significantly the fission reaction kicks loose a few neutrons that may
In this chapter we will review tfiose aspects of nuclear physics that are particu-
then go on to induce more fission reactions. Hence we can use the neutrons to
larly relevant to the study of fission chain reactions. It should be stressed that this
12 / INTRODUCTORY CONCEIT3 OF NUCLEAR POWER R E A a O R ANALYSIS THE NUCLEAR PHYSICS OF FISSION CHAIN REACTIONS / 13
presentation is not intended to be complete. Indeed we would anticipate that most naturally occurring nuclides include alpha decay, in which the nucleus emits a
engineering students will have had some exposure to modem atomic and nuclear helium nucleus :He; beta decay, which corresponds to the conversion of a neutron
physics in earlier courses. (Students who find most of the material in this chapter to in the nucleus into a proton, generally accompanied by the emission of an electron
be totally alien temtory would be well advised to consult one of the several and a neutrino; and gamma decay, the transition of a nucleus from one excited
excellent standard references on introductory nuclear physics1-" containing sub- state to a lower excited state with the accompanying emission of a photon.
stantially more thorough discussions of these topics.) Unfortunately, however, most However other types of radioactive decay are possible in a nuclear reactor since
conventional treatments of these subjects do not place sufficient emphasis on the many unstable nuclides are produced in fission which do not occur in nature. For
study of nuclear reactions in general, or neutron-nuclear reactions in particular, for example, certain nuclei such as !zKr may decay by emitting a neutron. (We will
our purposes (although there are several notable exceptionss). later find that this particular type of decay process is extremely important for
We will begin with a brief introduction to spontaneous nuclear radioactive decay reactor operation.)
as an example of a nuclear reaction. We then consider nuclear collision reactions The fundamental law describing radioactive decay is based on the experimental
and introduce the concept of a nuclear cross section. Here we will devote particular observation that the probability that a nucleus will decay in a given time-interval is
attention to a qualitative discussion of cross sections characterizing neutron- essentially constant, independent of the age of the nucleus or its environment,
nucleus reactions. Our final topic will be that of the nuclear fission reaction itself dependent only on the type of the nucleus itself. Hence the time rate of change of
and the radiation emanating from such reactions. the number of original nuclei of a given type must be proportional to the number
of nuclei present at that time. Let us call the proportionality constant A. Then if
N(t) is the number of original nuclei left at time I, we find
I. NUCLEAR REACTIONS
dN -AN (t).
dt
There are essentially two types of nuclear reactions of importance in the
study of nuclear reactors: (a) spontaneous disintegrations of nuclei and (b) reac- Here X is referred to as the radioactioe decay constant characteristic of the nucleus
tions resulting from the collision between nuclei and/or nuclear particles. An and has units of inverse time. If we initially have No nuclei present, then at any
example of the first type of reaction would be the radioactive decay of fission later time t the number of nuclei present will be given by an exponential law:
products, since these are frequently unstable. Such disintegration reactions depend
only on the properties of an individual nucleus. The neutron-nucleus collision
events involved in the fission chain reaction are an example of the second type of
reaction. These collision reactions depend not only on the properties of the The rate at which nuclei are decaying is given by
colliding particles, for example, the neutron and the nucleus, but also the relative
velocity with which they strike one another.
Before diving off into a discussion of nuclear reactions, let us first introduce
some notation. We will denote the number of protons in an atomic nucleus by Z From this time behavior, it is apparent that the probability that a given nucleus will
(the atomic number), the number of neutrons by N, and the total number of decay in a time interval I to t + dt is just
nucleons (protons plus neutrons) by A (the mass number). A specific nucleus will be
denoted by a symbol such as $X,where X is the chemical symbol for the atom of
interest. For example, iH, ':C, and 2U
: are notations for three such nuclei. We will
refer to various species of nuclei as nuclides. Nuclei characterized by the same Since radioactive decay is a statistical phenomenon, we cannot predict with any
atomic number Z but different mass numbers A are referred to as isotopes (e.g., certainty preci_selywhen a given nucleus will decay. However we can calculate the
2iiU, Z U , ,i:U, and :'U). Since the nucleus is a quantum mechanical system, it mean lifetime t of the nucleus before decay using our expression for p(t) from Eq.
may be found in any of a number of possible energy states. The general notation (2-4)
A,x refers to a nuclear ground state, while an asterisk is used to denote a nucleus in
an excited state, GX*.Long lived excited states of nuclei are referred to as nuclear
isomers or isomeric states and are denoted by a superscript m (e.g., in).
Hence on the average a given nucleus will decay after a time I/A.
A closely related quantity is the length of time necessary for half of the original
A. Radioactive Decay number of nuclei present to decay away. Such a time T,,, is referred to as
Certain nuclei are unstable in the sense that they may spontaneously undergo radioactive hag--life for the nucleus and can be calculated from its definition by
a transformation into a different nuclide, usually accompanied by the emission of noting
energetic particles. Such a spontaneous nuclear transformation is referred to as
radioactive decay. The three most common types of radioactive decay found in
-
- -
14 / INTRODUCTORY CONCEITS OF NUCLEAR POWER REACTOR ANALYSIS THE NUCLEAR PHYSICS OF FISSION CHAIN REACllONS / IS
concept here is the uncertainty or width r of the energy level characterizing the
excited state. This width is related to the mean lifetime of the state by the
Heisenberg uncertainty principle:
We can also write similar equations describing several nuclides, each of which
decays into another. Consider, by way of example, the radioactive decay chain:
Target Projectile
Then the appropriate equations describing the number of nuclides of each type
present are
As an example, the reaction
would be written as
where Rx(f) is the production term for the X-nuclide, and so on. Since this is just a The general class of such reactions would be simply denoted as ( n , y ) reactions.
system of linear first-order differential equations with constant coefficients, it can Nuclear reactions are generally accompanied by either the absorption or emis-
easily be solved using standard techniques, and hence we will defer further sion of energy. One can calculate the energy released by (or required for) a given
discussion to the problems at the end of the chapter. nuclear reaction by using the important result from the theory of relativity:
Very similar considerations also hold for the transition of nuclei between
different excited states. Such states represent the quantum levels available to the
nucleus. We can again characterize the probability that the nucleus will "decay" where c is the speed of light and m is the mass converted into energy in a reaction.
out of one excited state into a lower state by a decay constant ?, and once again The appropriate quantity to use for the variable m that appears in this formula is
also develop the concept of a mean lifetime for the excited state I. A useful related the mass difference between the interacting particles before and after the collision.
16 / INTRODUCFORY CONCEPTS OF NUCLEAR POWER REACTOR ANALYSIS THE NUCLEAR PHYSICS OF FISSION CHAIN REACTIONS / 17
For the reaction a(b,c)d we would calculate the reaction energy as ally by considering a beam of neutrons, all traveling with the same speed and
direction, which is incident normally upon and uniformly across the face of a target
of material. If the target is sufficiently thin (say, one atomic layer thick), then no
nuclei in the target will be shielded by other nuclei from the incident neutron beam
(see Figure 2-1). In this case we would expect that the rate of neutron-nuclear
If Q >O, then we say the reaction is exorhermic, which corresponds to a release of reactions in the target will be proportional to both the incident neutron beam
energy in the reaction. If Q<0, then the reaction is said to be endothermic, and intensity I (in units of number of neutrons/cm2-sec) and the number of target
energy must be supplied to the colliding nuclei in order to stimulate the reaction to atoms per unit area N, (#/cm2). If we call the constant of proportionality a, we
occur. Obviously, nuclear fission is an example of an exothermic reaction. can write the rate at which reactions occur per unit area on the target as
There are a wide variety of possible nuclear reactions. The reactions of most
interest in the analysis of a nuclear fission reactor involve interactions between
neutrons and nuclei and include
Nuclear fission (n, fission):
+ $X-P$X +$x + neutrons+ 200 MeV, (2- 15) We have indicated the units of each of these quantities since they imply that the
proportionality factor a must have the units of an area.
Radiative capture (n, y): If the incident neutrons and target nuclei could be visualized as classical
particles, a would quite naturally correspond to the cross sectional area presented
by each of the target nuclei to the beam. Hence a is known as the microscopic cross
section characterizing the probability of a neutron-nuclear reaction for the nucleus.
Scattering (n, n) or (n, n'): We might continue to think of a as the effective cross sectional area presented by
the nucleus to the beam of incident neutrons. Since the nuclear radius is roughly
An+$x+An+$~ [elastic scattering (n,n)] 10-l2 cm, the geometrical cross sectional area of the nucleus is roughly 10-" cm2.
Hence we might expect that nuclear cross sections are of the order of loFz4cm2. In
+An + ($X)* [inelastic scattering (n,nl)] (2- 17) fact microscopic cross sections are usually measured in units of this size called
-+An+$x + y [inelastic scattering (n, n')]. barns (b). However this geometrical interpretation of a nuclear cross section can
frequently be misleading since a can be much larger (or smaller) than the geometri-
cal cross section of the nucleus due to resonance effects which, in turn, are a
We have already discussed the nuclear fission reaction. In radiatiw capture the consequence of the quantum mechanical nature of the neutron and the nucleus.
incident neutron is absorbed by the target nucleus to form a new nuclide of mass
For example, the absorption cross section of :'Xe for slow neutrons is almost one
+
number A 1. As we will see later, this "compound" nucleus is formed in an million times larger than its geometrical cross section.
excited state. In a radiative capture reaction, it will eventually decay to its ground
state by emitting a high-energy photon, that is, a gamma ray. An alternative type of
capture reaction of some importance in reactor control is the (n,a) reaction which
occurs in '!B, for example.
-
The thud reaction of importance is scattering. In this reaction the neutron simply
scatters off of the nucleus (n,n), although in some cases, it may first combine with
the nucleus to form a compound nucleus for a short time before being reemitted
and will frequently leave the nucleus in an excited state from which it later decays 1 neutrondcm2~
sec
by gamma emission.
The importance of the fission reaction to nuclear reactor operation is obvious.
Both radiative capture and scattering are also extremely important since they
influence the neutron economy and hence the chain reaction. We will concentrate
specifically on neutron-nuclear reactions as we turn to a more quantitative treat-
ment of nuclear reactions of importance in fission chain reactions.
We can give a slightly more formal definition of the microscopic cross section by
rearranging Eq. (2-18) to write
In this sense, then, if the target has a total cross sectional area &, all of which is
uniformly exposed to the incident beam, then
a Probability per nucleus that a neutron
-P (2-20)
& in the beam will interact with it a, "8"
(elastic (inelastic
Thus far we have been discussing the concept of a nuclear cross section in a scattering) scattering)
rather abstract sense without actually specifying the type of reaction we have in
mind. Actually such cross sections can be used to characterize any type of nuclear
reaction. We can define a microscopic cross section for each type of neutron-
nuclear reaction and each type of nuclide. For example, the appropriate cross
sections characterizing the three types of reactions we discussed earlier, fissioh,
radiative capture, and scattering, are denoted by o,, a,, and a,, respectively. We can
also assign separate cross sections to characterize elastic scattering a, in which the FIGURE 2-2. Neutron rroar &on beburhy.
target nucleus remains in its ground state, and inelastic scattering a , in which the
target nucleus is left in an excited state. Since cross sections are related to
surface of a target. However it is certainly conceivable that such cross sections will
probabilities of various types of reactions, it is apparent that
vary, depending on the incident neutron speed (or energy) and direction. Indeed if
the mic~oswpiccross section for various incident neutron energies is measured, a
very strong energy dependence of the cross section is found. The dependence of
In a similar sense we can define the absorption cross section characterizing those neutron cross sections on the incident beam angle is usually much weaker and can
events in which a nucleus absorbs a neutron. There are a number of possible types almost always be ignored in nuclear reactor applications. We will return later to
of absorption reactions including fission, radiative capture, (n,a) reactions, and so consider in further detail the dependence of cross sections on incident neutron
on. (Actually one could argue that fission is not really an absorption reaction since energy. First, however, it is useful to develop a quantity closely related to the
several neutrons are created in the fission reaction. It has become customary, microscopic cross section, that of the macroscopic cross section.
however, to treat fission as an absorption event and then add back in the fission
neutrons released in the reaction at another point, as we will see later.) Finally, we 2. MACROSCOPIC CROSS SECTIONS
can introduce the concept of the total cross section a, characterizing the probability Thus far we have considered a beam of neutrons incident upon a very thin
that any type of neutron-nuclear reaction will occur. Obviously target. This was done to insure that each nucleus in the target would be exposed to
the same beam intensity. If the target were thicker, the nuclei deeper within the
target would tend to be shielded from the incident beam by the nuclei nearer the
surface since interactions remove neutrons from the beam. To account for such
A schematic diagram9 of the heirarchy of cross sections along with their conven- finite thickness effects, let us now consider a neutron beam incident upon the
tional notation is shown in Figure 2-2. Notice that in general one would define the
surface of a target of arbitrary thickness as indicated schematically in Figure 2-3.
absorption cross section to characterize any event other than scattering
We will derive an equation for the "virgin" beam intensity I (x) at any point x in
a1 = 0, - a,. the target. By virgin beam we are referring to that portion of the neutrons in the
beam that have not interacted with target nuclei. Consider a differential thickness
In a similar fashion, one occasionally defines a nonelastic cross section as any event of target between x and x + dx. Then since a3 is infinitesimally thin, we know that
other than elastic scattering the results from our study of thin targets can be used to calculate the rate at which
neutrons suffer interactions in d x per cm2. If we recognize that the number of
target nuclei per cm2 in d x is given by dN,= N d x , where N is the number density
Thus far we have defined the concept of a microscopic cross section by of nuclei in the target, then the total reaction rate per unit area in d x is just
c0nsidering.a beam of neutrons of identical speeds incident normally upon the
20 / INTRODUCTORY CONCEPTS OF NUCLEAR POWER REACTOR ANALYSIS THE NUCLEAR PHYSICS OF FISSION CHAIN REACTIONS / 21
Hence, it is natural to interpret Z, as the probability per unit path length traveled
that the neutron will undergo a reaction with a nucleus in the sample. In this sense
then
With this interaction probability, we can calculate the average distance a neutron
FIGURE 2-3. A monaac~@lcmm barn indent wrmrlly on I thick target.
h travels before interacting with a nucleus in the sample
Notice that, consistent with our prescription that any type of interaction will
deflower an incident neutron, we have utilized the total microscopic cross section 0,
in computing dR.
We can now equate this reaction rate to the decrease in beam intensity between
x and x + d x It is customary to refer to this distance as the neutron mean free path since it
essentially measures the average distance a neutron is likely to stream freely before
colliding with a nucleus.
The reader has probably noticed the similarity of this analysis to our earlier
Dividing by dx we find a differential equation for the beam intensity I ( x ) treatment of radioactive decay. The spatial attentuation of a neutron beam passing
through a sample of material and the temporal decay of a sample of radioactive
nuclei are similar types of statistical phenomenon in which the probability of an
event occurring that removes a neutron or nucleus from the original sample
If we solve this equation subject to an incident beam intensity of I, at x=O, we depends only on the number of neutrons or nuclei present at the position or time of
find an exponential attenuation of the incident beam of the form interest. It should be stressed that both the mean free path and the mean lifetime
for decay are very much average quantities. There will be statistical fluctuations
about these mean values.
If we recall that Z, is the probability per unit path length that a neutron will
The product of the atomic number density N and the microscopic cross section undergo a reaction, while the neutron speed o is the distance traveled by the
o, that appears in the exponential term arises so frequently in nuclear reactor neutron in a unit time, then evidently
studies that it has become customary to denote it by a special symbol: [cm-I]= [sec-'I= Frequency with which
reactions occur. (2-28)
This quantity is usually referred to as the collision frequency for the neutron in the
One refers to Z, as the total macroscopic cross section characterizing the target sample. Its reciprocal, [oZ,]-', is therefore interpretable as the mean time between
material. The term "macroscopic" arises from the recognition that Z, characterizes neutron reactions.
the probability of neutron interaction in a macroscopic chunk of material (the Thus far our discussion has been restricted to total macroscopic cross sections
target), whereas the microscopic cross section characterizes the probability of that characterize the probability that a neutron will undergo any type of reaction.
interaction with only a single nucleus. We can generalize this concept by formally defining the macroscopic cross section
It should be noted that X, is not really a "cross section" at all, however, since its for any specific reaction as just the microscopic cross section for the reaction of
units are inverse length. A more appropriate interpretation can be achieved by interest multiplied by the number density N characterizing the material of interest.
reexamining Eq. (2-22) and noting that the fractional change in beam intensity For example, the macroscopic fission cross section would be defined as
occurring over a distance dx is just given by
Notice also that are compiled, evaluated, and then organized into data sets to be used by nuclear
engineers. We will return in a later section to discuss in further detail such nuclear
X, = Z,+ Z,. cross section data sets. For convenience, however, we have included in Appendix A
a table of some of the more important cross sections characteristic of "thermal"
It should be stressed that while one can formally define such macroscopic cross neutrons-that is, of neutrons whose energies are comparable to the thermal energy
sections for specific reactions, our earlier discussion of neutron penetration into a of atoms in a reactor core at room temperature, E=0.025 eV-which serves to
thick target applies only to the total macroscopic cross section X,. We could not illustrate typical orders of magnitudes of these quantities.
extend this discussion, for example, to the calculation of the probability of neutron
penetration to a depth x prior to absorption by merely replacing X, in Eq. (2-24) by EXAMPLE: As a specific illustration, let us calculate the mean free path for a
Z,, since it may be possible for the neutron to undergo a number of scattering thermal neutron in graphite. According to the table in Appendix A, carbon is
reactions before finally suffering an absorption reaction. We can calculate these nearly a pure scatterer with microscopic cross sections 0 ~ ~ 4b. and
8 a,=4.0x loT3
specific reaction probabilities only after a more complete consideration of neutron b. If we assume a mass density of 1.60 g/cm3, we can calculate the atomic
transport in materials (Chapters 4 and 5). number density in graphite as N = .0803 X l d 4 ~ m - ~Consequently
. the macro-
The concept of a macroswpic cross section can also be generalized to homo- scopic scattering and absorption cross sections are
genwus mixtures of different nuclides. For example, if we have a homogeneous
mixture of three different species of nuclide, X, Y, and Z, with respective atomic
number densities N,, N,, and N,, then the total macroscopic cross section
characterizing the mixture is given by and the total cross section is
where af is the microscopic total cross section for nuclide X,and so on. It should The mean free path is therefore
be noted that such a prescription for determining the macroscopic cross section for
a mixture arises quite naturally from our interpretation of such cross sections as
probabilities of reactions.
As we mentioned earlier, all neutron-nuclear reaction cross sections (fission,
radiative capture, scattering, etc.) depend to some degree on the energy of the Notice how small the absorption cross section in graphite is compared to its
incident neutron. If we denote the neutron energy by E, we acknowledge this scattering cross section. Indeed since Z,/Z,= 1.2 x I d , one can infer that a thermal
dependence by including a functional dependence on E in the microscopic cross neutron in graphite will make some 1200 scattering collisions on the average before
section o(E) and hence by inference also in the macroscopic cross section X(E). being absorbed. This very low absorption cross section makes graphite an ideal
However the macroswpic cross section can depend on additional variables as material for nuclear reactor applications.
well. For example, suppose that the target material does not have a uniform
composition. Then the number density N will depend on the position r in the C. Characteristics of Neutron-Nuclear Cross Sections
sample, and hence the macroscopic cross sections themselves will be space-
Before considering in detail the various types of neutron-nuclear reactions
dependent. In a similar manner, the number densities might depend on time-
significant in nuclear reactor analysis, it is useful to give a brief discussion of some
suppose, for example, that the nuclide of interest was unstable such that its number of the relevant physics underlying the behavior of these cross sections. There are
density was decaying as a function of time. Therefore in the most general case we two aspects involved in the analysis of neutron cross cections:
would write (a) the kinematics of two-particle collisions and
(b) the dynamics of nuclear reactions.
The kinematics of two-body collisions, that is, the application of the laws of
to indicate the explicit dependence of the macroscopic cross section on neutron conservation of momentum and energy to such collisions, should be very familiar
energy E, position r, and time r . to the reader from introductory courses in mechanics or modem physics. However
In summary then, nuclear cross sections can be used to characterize the probabil- consistent with our attempt to make this presentation as self-contained as possible,
ity of various types of neutron-nuclear reactions occurring. They obviously will be and being well aware of the short half-life of this type of information in the fast
a very basic ingredient in any study of fission chain reactions. The determination core memory of most students, we will review such kinematic calculations in
of such cross sections is the task of the nuclear physicist and involves both Section 2-I-D.
experimental measurement and theoretical calculations. The enormous amount of The dynamics of nuclear reactions is concerned with the fundamental physical
cross section information required for nuclear reactor analysis is gathered by mechanisms involved in such collision events. The two mechanisms of most interest
numerous nuclear research centers throughout the world. These cross section data in nuclear reactor applications are those of potential scattering, in which the
--
24 / INTRODUCTORY C O N C E W OF NUCLEAR POWER REA(TOR ANALYSIS THE NUCLEAR PHYSICS OF FISSION CHAIN REACTIONS / 25
neutron merely bounces off of the force field of the nucleus without actually compound nucleus is much higher due to the additional binding energy of the
penetrating the nuclear surface, and compound nucleus formation, in which the added neutron E,. If E,+ Eb is very close to a nuclear energy level of the
incident neutron is actually absorbed by the nucleus to form a new nucleus of mass compound nucleus *:'x, one expects that the probability for compound nucleus
number A + 1 which then decays by emitting a gamma, a neutron, or perhaps formation will be much larger than if Ec+ Eb does not "match" this energy level.
fissioning. Hence we expect that the cross sections for such compound nuclear reactions will
We will give a brief discussion of both of these interaction mechanisms, and then exhibit sharp peaks or resonances at those neutron energies E for which this energy
turn our attention to a specific survey of the various types of cross section behavior matching occurs. By way of illustration, we have shown the resonance structure in
encountered in the more common materials utilized in nuclear reactors. the low-energy cross section behavior of a number of nuclei in Figure 2-4.
1. MECHANISMS OF NEUTRON-NUCLEAR INTERACTION
The simplest type of nuclear reaction occumng in a nuclear reactor is
potential scattering, in which the neutron scatters elastically off the nuclear poten-
tial without ever penetrating the nucleus itself. This type of collision event is very
similar to that which would occur between two hard spheres (e.g., billiard balls),
and the cross section for such a reaction is essentially just the geometrical cross
section of the nucleus. Potential scattering cross sections are characterized by a
rather flat energy dependence from about 1 eV up to the MeV range.
Another common type of reaction encountered in a nuclear reactor is that in
which the incident neutron is first absorbed by the nucleus $X to create a new
compound nucleus A ~ I xThis
. compound nucleus subsequently decays by emitting
an energetic particle. Compound nucleus formation occurs in many neutron-
nuclear reactions of interest to the reactor engineer, including fission, radiative
capture, and certain types of scattering.
That such a mechanism must be involved in these reactions can be inferred from
the relatively long times (at least on a nuclear time scale) for such events to occur.
Whereas it would take a slow neutron (traveling at Id cm/sec) some 10-l7 sec to
cross the nucleus, neutron-nuclear reactions such as fission occur on a time scale
of some lo-" sec--Or some 1000 transit times. Hence the incident neutron must
first be absorbed by the original nucleus and rattle around a bit, distributing both
its kinetic energy and the additional binding energy supplied by the added neutron
to the other nucleons in the nucleus, before the compound nucleus finally decays.
The long lifetime of the compound nucleus implies that the disintegration process
is essentially independent of the original mode of formation; that is, the compound
nucleus lasts long enough to "forget" most of the characteristics of the incident
neutron (such as which direction it came from).
The formation of a compound nucleus actually corresponds to a so-called
resonance reaction, in which the incident neutron energy matches one of the energy
levels in the compound nucleus. To be more specific, consider a neutron incident
upon a nuclide $X:
L I I I J
As we will see later when we develop the topic of collision kinematics, the energy 0.001 0.01 0.1 1.0 10.0
available for such a reaction is the center of mass (CM) energy E c = ( M / m + M ) E , Energy lev)
where m is the neutron mass, M is the nuclear mass, and E is the neutron kinetic
energy in the laboratory system. The actual energy of the excited level of the FIGURE 2-4. Low-energy cross W l o n behnvlor of several important nucllde~.'~
26 / INTRODUCTORY CONCEITS OF NUCLEAR POWER REACTOR ANALYSIS THE NUCLEAR PHYSICS OF FISSION CHAIN REACTIONS / 27
The second stage of a compound nucleus process is the decay of the compound
nucleus. This decay may occur in a variety of ways, as indicated schematically
below:
Since we have argued that the final disintegration mode is essentially independent
of the neutron absorption process which creates the compound nucleus, we might
expect that the cross section energy dependence for compound nucleus reactions
will exhibit certain similarities. These similarities will become apparent as we
consider in more detail the specific neutron-nuclear reactions of importance in
nuclear fission reactors.
Neutron in
t
I
7 Emission
Eo kc
FIGURE 2-6. A dngk-kvd capture resonance.
Neutron in
Neutron
mc' + M,C'
Neutron energy, eV
Since this wavelength decreases with increasing neutron energy, it is apparent that TABLE 2 1 A Summary d Gloss Sectbn Behnvlor for Nuekl of Varlous M.ss Numberno
for sufficiently high energies, the probability of neutron interaction with the
nucleus will similarly decrease. Such very high energy behavior has little relevance
Slow Epithermal Fast
to nuclear reactor analysis, since neutrons in fission chain reactions rarely exceed neutrons neutrons neutrons
IOMeV in energy. E<leV leV<E<.IMeV .lMeV< E <20MeV
The irregular, jagged behavior at low energies in region @ is also a wavelength I
effect For sufficiently small energies, the wavelength becomes comparable to the
interatomic spacing, and the neutron interacts not with a single nucleus but rather
with an aggregate of nuclei. If the material has a regular structure (such as the Separated resonances
Light
crystalline structure of graphite), the neutron will be diffracted, just as X-rays are nuclei Potential scattering
diffracted when passing through a crystal. This is accompanied by a sensitive A <2S t
energy dependence as the neutron wavelength becomes comparable to multiples of Resonance scattering, (n.2n). (n,p)
the spacing between various crystal lattice planes. For sufficiently small energies, Separated Overlapping Continuum
the wavelength becomes so large that diffraction becomes impossible, and the cross resonances
section becomes smoothly varying again. Intermediate
There are two other important effects which influence the neutron cross section nuclei Resonance scattering, radiative capture
behavior at very low energies (regions @ and 0). If the neutron energy is less
25<A<80
Potential scattering Inelastic scattering
than the chemical binding energy of atoms in the sample (- l eV), the neutron will
no longer be interacting with a free nucleus, but rather with an aggregate of bound <
nuclei. It can interact and excite the internal modes of the sample, such as crystal
lattice vibrations or molecular rotations. For these low energies, the thermal motion Separated Overlapping Continuum
of the nuclei also becomes very important. If we recall that these thermal motions mnance* resonances
Heavy
are essentially characterized by the sample temperature T and that the atoms in nuclei
thermal equilibrium at this temperature are characterized by a thermal energy f kT, A>80
where k is the Boltzmann constant, k = 8.6173 x 10-'eVl°K, then for neutron
energies comparable to this thermal energy (at room temperature, E = 0.025 eV), the
motion of the nuclei must be considered. We will return in the next section to
discuss the modifications this requires in the neutron cross section. Suffice it to say
t- Inelastic scattering, (n, 2n)
at this point that such considerations imply that for small neutron energies, the
cross section will behave as 1/ E 'IZ. Basic nuclear data appear in a wide variety of forms. For example, raw cross
Similar features are found in most neutron cross sections of interest in reactor section data are usually provided by a variety of experiments, each char-
analysis, although the particular energy regions in which such behavior arises will acterized by a different degree of accuracy. Multiple sets of experimental data may
vary from nuclide to nuclide. For example, we have summarized the characteristics exist giving different values for the same cross sections. Data may be provided by
of nuclei of different mass numbers in Table 2-l.1° different theoretical calculations of varying accuracy. There are also numerous
We will return later from time to time to discuss more specific features of gaps where no cross section measurement or theory is available or applicable. The
neutron cross section behavior which are of particular significance for nuclear enormous volume of such varied nuclear data would overwhelm the nuclear
reactor applications. engineer in his efforts to extract those cross sections of relevance to his particular
needs.
3. NUCLEAR DATA SETS Hence a number of years ago it was decided to consolidate and standardize all of
The key ingredient in all reactor calculations is a knowledge of the various the cross section information into one data set. To this end the Evaluated Nuclear
relevant neutron-nuclear cross sections. The complicated dependence of such cross Data File (ENDF)I4 was established to consolidate, organize, and present these
sections on neutron energy and angle of incidence, combined with the large data in a form convenient for nuclear applications. The ENDF system contains
number of isotopes involved in nuclear reactor analysis implies that neutron cross both neutron and photon cross section data along with data-processing computer
section data can be quite massive. Such data have been accumulated over the past programs which can manipulate the data into the most convenient form for the
few decades by both experimental measurements and theoretical calculations. The user. These data are stored in three computer library systems:
tabulation of these data has grown from a single volume, BNL325" (the so-called (1) CSISRS (Cross Section Information Storage and Retrieval System): This
"barn book") in the 1950s to a six-volume set'' in the 1960s to the point where data set contains essentially unevaluated raw data from experimental
nuclear data are most conveniently kept on magnetic tape or in computer memory. measurements.
(Although it should be mentioned that a new third edition" of BNL-325 has (2) ENDF/A: This data set contains both complete and incomplete sets of
recently been issued.) nuclear data as soon as they become available. For each isotope there
36 / INTRODUCTORY CONCEPTS OF NUCLEAR POWER REACTOR ANALYSIS THE NUCLEAR PHYSICS OF FISSION CHAIN REACIlONS / 37
space coordinates (0,cp) (see Figure 2-10). Notice that to describe- the incident collision changes the original neutron energy from E to E' in dE'. It is important to
neutron velocity we now specify both its energy E and its direction 51. notice that this differential scattering cross section is a "distribution" in the sense
At this point it is convenient to consider how one integrates over these variables that it is associated with a certain range of final energies, E' to E'+ dE'. Hence its
(as we shall have cause to do later in this section). Suppose we wished to integrate dimensions are cm2/eV.
over all possible neutron velocities. This integration could then be performed in There is a very simple relationship between the differential scattering cross
either Cartesian or spherical velocity coordinates: section o,(E+E1) and our earlier definition of the microscopic scattering cross
section o,(E). If we recognize that the latter quantity is just related to the
probability that a neutron of energy E will suffer a scattering collision, regardless
of the final energy E' to which it is scattered, then it is apparent that o,(E) is just
the integral of the differential scattering cross section o,(E+Ef) over all final
yowever we have defined the unit vector in the direction of the velocity vector v as energies E '
a. Hence we can identify the angular portion of the integration in Eq. (2-42) as just
the integration over this direction:
Such differential cross sections are quite important in nuclear reactor analysis
angle, p,,~cosO, which can bepo?veniently expressed as the dot product between since they determine the manner in which neutrons move about in a reactor core,
the unit direction vectors IL~-i2.62'. (Again, the dependence on azimuthal angle 9 as well as the rate at which they leak out of the reactor. The measurement or
does not arise for materials in which the nuclei have a random orientation.) One calculation cf such cross sections can become quite involved, and the amount of
occasionally denotes this functional dependence by writing data necessary to adequately represent differential cross sections is usually rather
voluminous. Such data are contained in evaluated cross section files such as
ENDF/B, as well as in cross section compilations such as B N L - ~ o o . ' ~
Although the calculation of such differential scattering cross sections is usually
We will continue to use the somewhat more formal notation by writing a , ( b + b l ) , formidable, there is one instance of considerable importance to nuclear reactor
even though we know that this differential cross section usually depends only on analysis in which such cross sections can be calculated in a straightforward fashion
16. merely by using the laws of conservation of energy and momentum. This is the
Thus far we have developed the concept of differential scattering cross sections situation in which neutrons scatter elastically from stationary nuclei. To prepare
that characterize the probability of scattering from one energy to another or one the way for the calculation of such cross sections, let us first decompose the
direction to another. We can combine these concepts by defining a double dgjeren- d:fferential scattering cross section into two factors:
tial scattering cross section that characterizes :catteripg from an incident energy E,
direction 52 to a final energy E' in dE' and n' in dn'
If we recall our earlier definition of the scattering cross section, u,(E), then it
becomes apparent that we can identify
where we have introduced the reduced mass p~ mM/(m + M). Hence we find the
important relation between the energy in the CM and the LAB frames as
Target
Neutron COM nucleus
-(8-
v vCU
@
m M In particular it should be noted that the total energy in the CM system is always
less than that in the LAB system. The energy difference is taken up by the center of
v\;
mass motion itself.
LAB - Before LAB - After Using conservation of momentum and energy, it is easy for one to demonstrate
that the magnitudes of the CM velocities do not change in the collision:
Scattaring
angle - CM
syltern
-Me--
only their velocity vectors are rotated through the CM scattering angle 0,. This fact
allows one to relate the scattering angles in the LAB and CM frames. Consider the
vector diagram in Figure 2-13 illustrating the velocities and scattering angles in
these two frames.
If we note from this diagram that
CM - Before CM - After vLsinBL= v&sinBc,
FIGURE 2-12. Detinitioml of d l l s l o n coordlmtes in LAB and center-of-m8ss (CM) systems. 0; cos BL = V +
C ~ 0;: cos Bc,
then we can relate the scattering angles in the CM and LAB frames by
to the neutron and upper-case notation to the nucleus. The subscripts L and C refer
to LAB or CM frames, respectively. vk sin Bc sin 0,
tan 9, =
The velocity of the CM frame is defined by
o, + 0;: cos Bc = -1 +cosBc
A
This relationship is particularly useful since cross sections are usually calculated
in the CM frame, but are measured and used in the LAB frame. If we denote the
where we have assumed that the initial nucleus velocity V, is zero and noted that differential scattering cross sections characterizing scattering through angles 9, and
the nucleus-neutron mass ratio M/m is essentiallyjust the nuclear mass number A. 8, in the LAB and CM frames, a,(@,) and oC,(Bc) respectively, we can use
If we note that the neutron and nucleus velocities in the CM frame are given by
oL(BL)sinBLdBL=aCM(BC)
sin flcdBc (2-62)
then it is apparent that the total momentum in the CM frame is zero, as it must be.
We can relate the total kinetic energies in the LAB and CM frames by
computing
0
LAB: &=fmv:+
CM: Ec= ;mu;+ ~ M V ; =f
FIGURE 2-13. Relation between the scattering an* In the LAB m d CM frames.
42 / INTRODUCTORY CONCEPTS OF NUCLEAR POWER REAmOR ANALYSIS THE NUCLEAR PHYSICS OF FISSION CHAIN REACTIONS / 43
to relate the LAB and CM differential scattering cross sections by If we recall that there is a one-to-one relationship between the neutron energy
transfer and the scattering angle, we can infer that there must similarly be a
relationship between the probability of the neutron experiencing a given energy
transfer Ei+E, and the probability that it will be scattered through a given
scattering angle 8,. Yet the probability of scattering through an angle 8, into doc
about Oc is just given by
Returning to our vector diagram in Figure 2-13 and using the law of cosines, we
can find
Notice in particular that the probability of scattering from an energy Ei to a final lattice position by the collision, but will possess sufficient recoil energy to dislocate
energy E, is independent of the final energy E,. other nuclei in the lattice, leading to significant radiation damage to the material.
Using this probability distribution, we can compute the average final energy of a This is an extremely important process in materials exposed to the high radiation
neutron suffering an elastic scattering collision as environment of a nuclear reactor core and must be taken into account in nuclear
reactor design. We will return to consider it in more detail in Chapter 11.
velocities A second example of interest is that in which the cross section a(lv-VI) is a very
slowly varying function of relative speed. Since the nuclear velocity distribution is
%(v)d3vZ Number of target atoms/cm3 with velocities rather sharply peaked around V = V,,, we can approximate Eq. (2-80) by neglecting
(2-78)
V in the relative speed in the integrand for neutron speeds v> V,,. In this case, we
V in d3V about V.
again find that
Then the probability of a neutron interaction per neutron per second with a
nucleus of any velocity is obtained by averaging Eq. (2-77) over %(V)
such that the true cross section and the observed cross section are again the same,
and the observed cross section is temperature-independent. This behavior is
frequently exhibited by the scattering cross section for many reactor materials.
Notice that we have essentially used this expression to define an averaged cross It is also of interest to examine the behavior of the averaged cross section for
section B(o) depending only on the neutron speed o. That is, we can write the small neutron speeds o& V,. Then we can replace Iv-V( by V in the integrand to
appropriately averaged cross section characterizing neutrons moving with a speed o find
through a sample of nuclei with a velocity distribution %(V) as
Since the integral is now just a constant, independent of the neutron speed o, we
It is important to remember that this is the cross section that would be measured in find that the behavior of S(o, T) as the neutron speed o becomes very small is just
an experiment (such as the transmission experiment we described earlier). No
experiment looks 'directly at the true neutron-nuclear reaction cross section, but
rather measures an average of a((v - VI) over the distribution of nuclear velocities.
The same is true for applications of such cross sections to nuclear reactor analysis.
This explains the low-energy behavior we observed earlier in the total cross section
Since any reactor is at a finite temperature, the nuclei comprising its core will be in
thermal motion, and hence one must take care to always use cross sections that for graphite. Of course this result can be explained physically by simply recogniz-
ing that for small neutron velocities, the neutron appears to be essentially
have been appropriately averaged over these nuclear velocities. Since this distribu-
stationary to the more rapidly moving nuclei. Hence the collision rate ceases to
tion function depends on the temperature characterizing the material of interest,
the "thermally averaged" cross sections will similarly depend on temperature. depend on the neutron speed and depends only on the nuclear speed distribution
It is useful to consider the application of this result to two particularly simple (i.e., the temperature of the scattering material).
examples of cross section behavior. We mentioned earlier that many nuclear cross Such averages over the thermal velocity distribution of the nuclei in a material
sections behave essentially as I/o (e.g., below a capture resonance). If nuclear must always be performed in measuring or utilizing cross sections characterizing
motion is to be included, we would express such behavior as thermal neutrons. However such an average is also extremely important in deter-
mining the correct effective cross sections to use when describing resonance
Y behavior.
a((v- V))= - Life could become exceedingly complex indeed if a detailed estimate of the
(v-VI '
nuclear velocity distribution function %(V) were required, since this depends on
If we substitute this form into Eq. (2-80) and note that the nuclear velocity the complicated microscopic dynamics of atoms in the reactor (e.g., atomic
distribution function is normalized such that vibrations in crystalline lattices or atomic motions in liquids). Fortunately it is
sufficient for most purposes to represent the nuclear velocity distribution by the
Maxwell-Boltzmann distribution characterizing an ideal gas in thermal equilibrium
at a temperature T :
then we find the averaged cross section becomes just
Hence in this instance the observed cross section depends on neutron speed o in
exactly the same way as the true cross section depends on relative speeds.
Furthermore the observed cross section is independent of temperature. We will make use of this particular average in the next section.
48 / INTRODUCTORY CONCEPTS OF NUCLEAR POWER REACTOR ANALYSIS THE NUCLEAR PHYSICS OF FISSION CHAIN REACI'IONS / 49
(b) THE DOPPLER EFFECT ON CROSS SECTION RESONANCE BEHAVIOR cast Eq. (2-91) into a slightly different form by defining the variables
We have seen that the actual or true cross section depends on the relative
speed between the neutron and the target nucleus. However since the nuclei x r 2 ( E - E o ) / r and J =I./TD
themselves are in thermal motion, this relative speed may be either greater or less
than the neutron speed. This difference in relative speeds gives rise to a "Doppler where TD is the so-called Doppler width of the resonance
shift" effect in resonance cross section behavior.
We can take account of this effect by merely substituting the Breit-Wigner
resonance cross section formulas into our expressions for the thermally averaged
cross sections in Eqs. (2-80) or (2-88). For the purposes of this calculation it is
usually adequate to assume that the nuclear velocities are described by a Maxwell- Then one can write
Boltzmann distribution M (V).
Now recall that the single-level Breit-Wigner formula for a capture resonance
gives the cross section
where
in terms of the center of mass energy Ec Our task is to express E, in terms of the
nuclear velocity vector V and then perform the integration over nuclear velocities In practice, the * ( J , x ) would be calculated for each value of J and x of interest
indicated in the averaging formula Eq.(2-88). There is really nothing particularly using straightforward numerical integration techniques. To better understand the
complicated about this task, except for a bit of vector algebra and the fact that the implications of such calculations, we have sketched the thermally averaged capture
integral that arises cannot be explicitly performed (rather, it must be performed cross section as determined by Eq. (2-94)in Figure 2-14.
numerically or tabulated). In particular, we have sketched the dependence of this cross section on energy
Since this manipulation is not particularly enlightening, we will only outline the for several different temperatures T . It should first be noted that as the temperature
major steps here and refer the interested reader to more exhaustive treatments that T increases, the resonance broadens, while its peak magnitude decreases. For this
exist in numerous places throughout the We begin by noting that for reason, one frequently refers to resonance cross sections that have been averaged
the case of a Maxwell-Boltzmann distribution of target nuclei, one can partially over the distribution of nuclear velocities as "Doppler-broadened" cross sections.
performm the integration over nuclear velocity V to rewrite Eq. (2-88) as It should be stressed that we still have not introduced any additional assump-
tions or approximations into this derivation of the Doppler-broadened resonance
cross section form Eq. (2-94). And in several modern computer codes, Doppler-
(2-90)
where v, = lv - Vl while v , f( k ~ / m ) ' /If~we
. substitute the Breit-Wigner formula
(2-30) for or(v,) into this Integral, we find an exact expression for the averaged
cross section
(a) Neglect the second exponential in the integrand of \k({,x). 0.05 0.04309 0.04308 0.04306 0.04298 0.04267 0.04216 0.04145 0.04055 0.03380 0.01639
@) Replace 0'- u; by (u"- u19/2u'-which is equivalent to approximating 0.10 0.08384 0.08379 0.08364 0.08305 0.08013 0.07700 0.07208 0.06623 0.03291 0.00262
0.15 0.12239 0.12223 0.12176 0.1 1989 0.1.1268 0.10165 0.08805 0.07328 0.01695 0.00080
0.20 0.15889 0.15854 0.15748 0.15331 0.13777 0.11540 0.09027 0.06614 0.00713 0.00070
0.25 0.19347 0.19281 0.19086 0.18324 0.15584 0.11934 0.08277 0.05253 0.00394 0.00067
0.30 0.22624 0.22516 0.22197 0.20968 0.16729 0.11571 0.07042 0.03880 0.00314 0.00065
0.35 0.25731 0.25569 025091 0.23271 0.17288 0.10713 0.05724 0.02815 0.00289 0.00064
(c) Extend the lower limit of integration t o y = - m; 0.40 0.28679 0.28450 027776 0.25245 0.17359 0 . W 0.04566 0.02109 0.00277 0.00061
0.45 0.31477 0.31168 0.30261 0.26909 0.17052 0.08439 0.03670 0.01687 0.00270 0.00064
then one finds 0.50 0.34135 0.33733 0.32557 0.28286 0.16469 0.07346 0.03025 0.01446 0.00266 0.00063
The x Function
while
The Bethe-Placzek cross section Eq. (2-98) is used very frequently in reactor
calculations and for convenience we have tabulated #(J,x) in Table 2-2. The
approximations used to obtain the Bethe-Placzek form break down for high-
temperature target distributions and for low energy resonances. For example, a
comparison for the 0.296 eV fission resonance of U9Pu at T52000°C shows a
difference in broadening of 50% between the exact and approximate fonnula-
tions.=
However since it is still traditional to utilize the Bethe-Placzek form in analyzing t ~ D.. Beynon and I. S. Grant, NucL Sci Eng. 17, 547 (1963).
resonance behavior, we will include a brief analytical study of its behavior along
with that of the exact expression. First notice that for low temperatures T+O we
can see from Eqs. (2-92) and (2-93) that l+oo. Hence the integrand of both \k(J,x) At the other extreme, the high temperature limit T+oo implies that 1 4 0 . In this
and +(l,x) will vanish except in the neighborhood of y-x. We can therefore case we find
replace y by x in the denominator of the integrands to write
which is a Gaussian shape characterized by the Doppler width TD rather than the
"natural" line width r. Hence as the temperature increases, the resonance broadens
out from its natural width to eventually approach a width that depends on the
which is a comforting, although certainly expected, result. temperature as Tf.
52 / INTRODUCTORY CONCEPTS OF NUCLEAR POWER REACTOR ANALYSIS THE NUCLEAR PHYSICS OF FISSION CHAIN REACTIONS / S3
One other observation is important. If the BethePLaczek form is examined in (c) DIFFERENTIAL SCATTERING CROSS SECTIONS WITH UPSCATTERING
greater detail, it becomes apparent that regardless of the temperature, the area We will now discuss the modifications necessary in our study of the differen-
under the resonance remains constant. We can demonstrate this by simply integrat- tial scattering cross section characterizing elastic potential scattering when effects
ing over the resonance of nuclear motion must be taken into account. Recall that in our earlier discussion
we found that the neutron could not gain energy in an elastic collision with a
stationary nucleus. It can only lose energy in such a collision. It is customary to
refer to such a process as "downscattering," since the neutron will scatter down in
energy.
However since the contribution to the integral from the "wings" of the resonance We can make this more explicit by considering a particularly simple example in
(far-off resonance) are so small, we can approximately extend the range of which an incident neutron scatters elastically from a stationary hydrogen nucleus
integration to 1+- co so that we can explicitly perfonn the integral to find (a proton). Then the scattering probability distribution is just
I , E l < Ei
P (Ei+ E,) =
0, El> E,
a result which is temperature-independent. That is, the scattering probability is independent of the final energy E, and
It should be pointed out that such behavior is a consequence of the Bethc vanishes for E, > E, (corresponding to upscattering in energy).
Placzek approximation. More generally, the area under the Doppler-broadened Let us now consider the situation in which the neutron suffers elastic scattering
resonance given by the exact expression Eq. (2-94) will in fact change with collisions in a hydrogen gas at finite temperature T in which the nuclei are in
temperature." However for the temperature range and resonances of interest in motion with a Maxwell-Boltzmann velocity distribution M(V,T). It is necessary to
most reactor applications, this change is relatively small, and one can essentially repeat our earlier consideration of two-body kinematics to include the motion of
assume that the area under the resonance is relatively insensitive to temperature the target nucleus. One would then have to average the cross section characterizing
changes. This fact will prove of some importance in our later study of reactor such scattering over the nuclear velocity distribution M (V, T), just as we did in the
behavior, since it will imply that an increase in the temperature of the absorbing previous sections. Since these tasks are rather tedious, we will simply note that the
material will increase the rate at which neutrons are absorbed in the resonance results of such calculations2' are that for such a proton gas at temperature T, the
range of the cross section. differential scattering cross section is given by
Thus far we have only examined the effect of thermal nuclear motions on a
capture resonance. However we could of course have also substituted in our
expression for scattering cross sections in the vicinity of a resonance, Eq. (2-38),
into the expression for thermally averaged cross sections Eq. (2-88) and labored
through some algebra to again arrive at Doppler-broadened cross sections which
again involve integrals that cannot be performed explicitly. We will avoid the
agony of such manipulation and only state the result of such labors here:
where we have defined another tabulated function2' x({,x) (see Table 2-2) to
characterize the interference term:
We have plotted P(Ei-+E,) for several incident neutron energies Ei in Figure 2-15.
It should first be noted that unlike the situation in which the hydrogen nuclei were
initially at rest, the scattering probability now depends on the final energy E,.
We will later utilize these expressions to study the absorption of neutrons in Furthermore this probability is not zero for E l > Ei for a finite temperature gas,
nuclear reactors. Although we have consistently denoted the thermally averaged hence implying that it is possible for the neutron to gain energy in a scattering
cross sections with an over-bar, for instance, C,(E,T), we will usually omit this collision. Such "upscattering" events are significant for incident neutron energies
notation in future discussions and simply regard all resonance cross sections as up to about IOkT. Above this energy, the scattering probability begins to resemble
having been Doppler-broadened. that characterizing stationary nuclei (i.e., a T=O hydrogen gas).
9 / INTRODUCI'ORY CONCEPTS OF NUCLEAR POWER REACTOR ANALYSIS THE NUCLEAR PHYSICS OF FISSION CHAIN REACTIONS / 55
have shown the fission cross sections characterizing the principal fissile nuclides,
2 3 3 ~ ,2 3 5 ~ ,and uspup taken from ENDF/B-IV.
The indicated cross section behavior is very similar to that of radiative capture
cross sections. However this would be expected since we have seen that compound
nucleus formation via neutron absorption is essentially independent of the mode of
compound nucleus disintegration or decay, for example, via fission or gamma
emission. It is particularly important to note that the fission cross section is over
two orders of magnitude larger for low-energy or thermal neutrons than for
high-energy fast neutrons (above I keV): The thermal neutron fission cross sections
are indeed enormous for these fissile isotopes, ranging up to thousands of barns in
magnitude. Such behavior will prove of very considerable importance in our later
studies of nuclear reactors.
We have also indicated the fission cross sections characterizing the principal
fissionable nuclides of interest, '"Th, 2 3 8 ~ , and 14"pU (see Figure 2-18). This
cross section behavior is somewhat different than that characterizing fissile nuc-
lides since fissionable nuclides can only be fissioned by sufficiently high-energy
neutrons. This implies that their fission cross sections will have a threshold energy,
below which the cross section drops to zero. Even above this threshold energy
(roughly 1MeV), the fission cross sections are quite low, being less than two barns.
When a neutron is absorbed by a fissile isotope such as u5U, it may induce that
isotope to fission. Yet it is also possible that the compound nucleus formed by the
neutron absorption, u6U*, might simply decay to its ground state by gamma
emission. The relative balance between the probability of fission and radiative
capture is an extremely important factor in nuclear reactor applications. We
characterize this balance by the capture-to$ision ratio, defined by
This ratio depends not only on the isotope of interest, but as well on the incident Yet just as significant as the energy released in the fission reaction is the fact that
neutron energy E. It is plotted in Figure 2-19 for the three primary fissile nuclides. several neutrons are also produced in the reaction. These neutrons can be used to
It can be seen that most neutron absorption in such isotopes leads to fission events propagate a fission chain reaction. Most of these fission neutrons appear essentially
(with the exception of a small range of a > l for 23'U).It should be noticed that a instantaneously (within 10-14 sec) of the fission event. These neutrons are referred
decreases quite appreciably above 0.1 MeV. This latter fact will prove to be of to as prompt. However a very few neutrons (less than 1%) appear with an
considerable importance when we discuss the concept of a fast breeder reactor. appreciable time delay from the subsequent decay of radioactive fission products.
Although only a very small fraction of the fission neutrons are delayed, thebe
C. Nuclear Fission Reactions delayed neutrons are vital for the effective control of the fission chain reaction.
The total number of neutrons (both prompt and delayed) released in a fission
A typical nuclear fission reaction such as reaction will vary. However in most nuclear applications we only need concern
ourselves with the average number of neutrons released per fission, which we
An + 2~&J+(~~U*)+fission
reaction products denote by v. This quantity will depend on both the nuclear isotope involved and
the incident neutron energy, generally tending to increase with increasing neutron
spews out a variety of reaction products, including the fissioned nuclei or fission energy. We have shown v(E) as a function of energy for the principal fissile
products and several neutrons as well as numerous gammas, betas, and neutrinos. isotopes in Figure 2-20.
And of course it also releases a very considerable amount of energy. Indeed a The neutrons produced in the fission reaction emerge with a distribution of
glance at the binding energy per nucleon before and after the fission reaction from energies, with the average fission neutron energy being roughly 2MeV. As with
Figure 2-16 suggests that energy on the order of 2OOMeV will be released in each othei fission parameters, this distribution will depend on the nuclear isotope
fission reaction. involved, to a lesser degree on the incident neutron energy, and will differ for
The fission fragment nuclei produced by the fission reaction are both highly prompt and delayed neutrons. To characterize this variation in fission neutron
charged and highly energetic. They slow down via collisions with adjacent atoms, energy, it is convenient to define the fission neutron energy spectrum, or more
losing energy and charge (picking up electrons) in the process. This is in fact the simply, the fission spectrum, x(E), defined as
principal mechanism by which the fission energy eventually appears as heat
generated in the fuel material. However these fission products are usually quite x ( E ) dE = Average number of fission neutrons emitted with
unstable as well, being somewhat neutron-rich, and will subsequently decay,
usually via beta emission. The energy released in such radioactive decay reactions
+
energy E in E to E dE per fission neutron. (2-1 11)
can amount to as much as 4-5% of the total energy released in the fission reaction.
Since such-"decay heat" will appear with an appreciable time delay corresponding
to the half-lives of the various nuclei involved, it can lead to difficulties unless
properly anticipated in fission reactor design.
-
0-
t-
a
0
13 U-235
A U-233
a X PU-233
1 OATR I N CLOSiLY-SPRCIG
Z RESW2NCE REGION SIPFRESSEO 1
29
m-
m
.-5
,-4
W
t o
c .
+ - ,- d p a
,a?:i
3
k-
(L
a
.=.
. ,
U
4
.C
.%.
t'
I
D
.',*.
0 i I I I I I I I I I I I I I I
lo-' 10-1 100 10'1 loeq lot5 1c.l 0 5 10 15
NEUTRON FNERGY. EV E (MeV)
FIGURE 2-19. Vdatlon of a with energy lor "'U. '"U, and U 9 ~ u . FIGURE 2-20. Avenge neutroa number per b b n Y as a lunetloa d enelby2'
62 / INTRODUCTORY CONCEPTS OF NUCLEAR POWER REACTOR ANALYSIS THE NUCLEAR PHYSICS OF FISSION CHAIN REACTIONS / 63
Since the relative isotopic yield per fission will vary for different fuel isotopes,
A typical prompt neutron fission spectrum is shown in Figure 2-21. One can also the detailed characteristics of the precursor groups will similarly be isotope-
represent the fission spectrum by a simple empirical expression; for example, the dependent. To this end, let us define:
prompt neutron fission spectrum of 235U is given by
&=Decay constant (P-decay) of ith precursor group
&=Fraction of all fission neutrons (both prompt and delayed) emitted per
fission that appear from ith precursor group
This proves useful in performing simple estimates of fission chain reaction be- p =Zi/3, =Total fraction of fission neutrons which are delayed.
havior. However more elaborate studies would simply use the fission spectra
tabulated in a nuclear data set such as ENDF/B. In Table 2-3 we have listed the half-lives, relative yield fractions Pi/P, and total
Because of the importance of delayed fission neutrons to nuclear reactor control, delayed neutron yields v, = vb, for the precursor groups characterizing the principle
it is useful to introduce a few related concepts useful for their description. By way fissionable isotopes recommended by ENDF/B-IV~'. Although dependent on the
of example, consider a typical fission product decay scheme leading to the emission fuel isotope, these data do not depend sensitively on the incident neutron energy
of a delayed neutron as sketched in Figure 2-22. It should be noted in particular below about 4MeV and hence can be used for either thermal or fast reactor
that the decay sequence leading to the delayed neutron emission is first the beta analysis.
decay of 07Br to "Kr*, followed by the subsequent decay of 87Kr* to via The energy spectrum of delayed fission neutrons is considerably lower than that
neutron emission. The effective time delay of this process is controlled by the of prompt fission neutrons and again depends on both the delayed neutron group
beta-decay-in this case, the half-life is some 55 sec. We refer to the fission and fissioning isotope. We have given a rough composite delayed neutron fission
fragment whose beta-decay yields a daughter nucleus which subsequently decays spectrum in Figure 2-23 along with typical measured spectra for the 55-sec and
via delayed neutron emission as a delayed neutron precursor. Of course a very large 22-sec groups of 235U. More detailed spectrum data can be found in the review
number (at least 45) of different delayed neutron precursor isotopes will be article of
produced in a fission chain reaction. It has been customary (and found to be Actually there are additional processes that can contribute delayed neutrons to
adequate) in reactor analysis to group these precursors into six classes the chain reaction. Photoneutron reactions (y,n) are particularly important in
characterized by approximate half-lives of 55, 22, 6, 2, 0.5, and 0.2-sec, respectively. reactors containing appreciable amounts of deuterium or beryllium. The decay
Each precursor group will contain a number of different isotopes. For example, times of these processes are even longer than those characterizing delayed fission
while the 55-sec precursor group is due almost entirely to one precursor, 8 7 ~ rthere
, neutrons (ranging up to 125 m).25 However these photoneutrons can usually be
are at least two major contributors, "Br and "'I, to the 22-sec group. The accounted for in reactor analysis by simply including one or more additional
composition of the remaining groups are considerably more complex. groups of delayed neutrons such as those tabulated above. Furthermore one can
64 / INTRODUCTORY CONCEPTS OF NUCLEAR POWER REACTOR ANALYSIS THE NUCLEAR PHYSICS O F FISSION CHAIN REACTIONS / 65
"%:
Group
-
rd 0.0527 i0.0040 n / j
Ti (sec) Relative Yield
"'U:
Group
vd -0.0074_+0.0004
Ti (scc) Relative Yield
235U:
vd-0.01668f 0.00070 n/f "'U: rd=0.046020.0025 n / j
Group Ti(sec) Relative Yield Group T, (sec) Relative Yield
Nonfission reactions
due to neutron
capture 4 100 cm delayed
66 / INTRODUCTORY CONCEPTS OF NUCLEAR POWER REACTOR ANALYSIS THE NUCLEAR PHYSICS OF FISSION CHAIN REACTIONS / 67
The majority of the fissior~energy appears as the kinetic energy of the fission 97% of the recoverable fission energy is deposited directly in the fuel material. The
fragments and is deposited essentially at the point of fission in the nuclear fuel. remainder is deposited in the coolant or structural materials by neutrons and
Note, however, that some of the fission energy appears as kinetic energy of gamma radiation, with less than 1% typically being deposited in shielding due to
neutrons (3%) and gammas (4%) with relatively long ranges. This energy will be gamma radiation. Actually as we will see in Chapter 12, the energy deposited in
distributed over the core of the reactor and adjacent material such as shielding. In other regions of the reactor is usually reassigned to the fuel in order to simplify the
Figure 2-24 we have noted the types of emergent radiation." thermal analysis of the reactor.
Furthermore it should be noted that some 4% of the fission energy appears in the
form of heat generated by the decay of radioactive fission products. If the nuclear
reactor were to be suddenly shut down, this decay heat would continue to be TABLE t S Z 6 Eflectlve Energy Rekued In Flssib.
produced and would have to be removed; otherwise the reactor core temperature The effective energy released in and following fission of the principal fissile isotopes by thermal
would rise dramatically, causing fuel element melting and failure. The removal of neulrons are:
such decay heat is one of the most serious problems in reactor safety studies.
Notice also that a sizable amount of energy (as much as 20MeV per fission) may lI3U: 190.020.5 MeV/f
be liberated by the high-energy gammas produced in radiative capture ( n , y ) 235U: 192.9e0.5 MeV/f
reactions. 239F'~: 198.5 5 0.8 MeV/f
It is customary to use an effective energy release per fission in determining the U'Pu: 200.320.8 MeV/f
portion of the total energy of fission that can be recovered by a coolant and hence The effective energy released following fission of the major fissionable isotopes by 235Ufission spectrum
contributes to the thermal power output of the reactor. Although this energy will neutrons are:
vary somewhat with the type of reactor and the detailed core composition, it is
typically of the order of 192MeV. (A more detailed tabulation of useful energy
release per fission has been calculated using atomic mass data by Jamesz6 and is
given in Table 2-5.) Of this 192MeV, some 168MeV appears as fission fragment
energy, while 7MeV appears as beta energy. These short-range contributions
deposit their energy in the nuclear fuel. If we also take into account the energy
deposited in the fuel (-7%) due to fast neutrons and gammas, we find that some
These values include all contributions except from neutrinos and very long-lived fission products.
D. Fission Fuels
Our previous discussion has indicated that there are a number of possibilities
available for fueling a fission chain-reacting system. In particular, we have noted
that the principal nuclidest of concern in nuclear reactor applications are:
f 4 neutrons
@.
@ /'
CJ
Neutron
slowing
down
4
/
Radiative
uoture
-
Capture 7's
Fissile nuclides: 233U,23'U, 2 3 9 ~2 ~4 ,1
while those susceptible to fast neutron fission are:
~ ~
Decay P's be fashioned out of natural uranium with even this low concentration of 235Uif one
*A bit of conventional notation2' for such nuclides remains as debris from the secrecy of the
Decay+ atomic weapons program during World War 11. Two-digit code numbers are used to identify
each isotope where the first digit is the atomic number minus 90 and the second digit is the last
FIGURE 2-24. Typs of radiation produced la a fluion &In reactloo. digit of the mass number. Hence "'U is denoted as "25." 238Uas "28." 239Puas "49," and so on.
61) / INTRODUCTORY CONCEPTS OF NUCLEAR POWER REACTOR ANALYSIS THE NUCLEAR PHYSICS OF FISSION CHAIN REACTIONS / 69
is sufficiently clever, most present-day reactor types are fueled with uranium in
which the percentage of 235U has been increased or enriched above its natural
value. As we will see later, such uranium enrichment is an extremely complicated
and expensive process.
There is yet another way to obtain fissile isotopes, however. It is found that when
certain nuclides absorb neutrons, they then undergo a sequence of radioactive
disintegrations that eventually result in the formation of a fissile isotope. The two
most important examples of such neutron transmutation reactions are:
Isotopes that can be transmitted into fissile nuclides via neutron capture are
referred to asfertile. The fertile isotopes of most interest are 238Uand n 2 ~ hwhich
,
are in abundant supply throughout the world.
Yet where does one find the neutrons necessary for this process? In a nuclear
reactor. Indeed since most present-day reactors are fueled with low-enrichment
uranium that may contain as high as 98% '"U, such transmutation processes will
occur quite naturally as the fertile nuclei capture excess neutrons from the fission
chain reaction. The key parameter in such processes is the number of neutrons
produced in each fission reaction per neutron absorbed in the fuel nuclei. (Here we
must remember that not all neutron absorptions in the fuel lead to fission--some
result in radiative capture.) We will define
REFERENCES
2-6. Boron is a common material used to shield against thermal neutrons. Estimate the 2-20. Determine the fission-rate density necessary to produce a thermal power density of
thickness of boron required to attenuate an incident thermal neutron beam to 0.1% of 400 kW/liter (typical of a fast breeder reactor core). Assume that the principal fissile
its intensity. (Use the thermal cross section data in Appendix A,) isotope is n 9 ~ ~ .
2-7. Suppose we consider a beam of neutrons incident upon a thin target with an intensity 2-21. An indium foil is counted at 5:00 p.m. Tuesday and found to yield 346,573 CPM in a
of 10'' neutrons/cml-sec. Suppose further that the total cross section for the nuclei in counter with a 50% efficiency for the 54-min In-116m activity. What is the probability
this target is 4 b. Using this information, determine how long one would have to wait, that none of these radioactive In-116m nuclei will remain in the foil at 2:00 p.m.
on the average, for a given nucleus in the target to suffer a neutron interaction. Thursday, the same week? (Note: ln(1- x)- - x,x< 1) [Victims working this problem
2-8. A free neutron is unstable against beta decay with a half-life of 11.7 m. Determine the can thank Dr. Ronald Fleming .)
relative probability that a neutron will undergo betadecay before being absorbed in 2-22. Compute and plot the parameter q for uranium enriched in 235U as a function of its
an infinite medium. Estimate this probability for a thermal neutron in HlO. enrichment (atom percent 2 3 5 ~at) thermal neutron energies.
2-9. Determine the number of scattering collisions a thermal neutron will experience on the
average before being absorbed in H20, DlO, 13", and cadmium, respectively.
2-10. How many mean free paths thick must a shield be designed in order to attenuate an
incident neutron beam by a factor of 1000?
2-11. Using the data from BNL-325, compute the mean free paths of neutrons with the
following energies in the specified materials: (a) 14MeV neutrons in air, water, and
uranium (characteristic of thermonuclear fusion neutrons), (b) 1MeV neutrons in air,
water, and uranium (fast breeder reactor neutrons), and (c) 0.05 eV neutrons in air,
water, and uranium (thermal reactor neutrons).
2-12. Determine the kinetic energy at which the wavelength of a neutron is comparable to:
(a) the diameter of a nucleus, (b) an atomic diameter, (c) the interatomic spacing in
graphite, and (d) the diameter of a nuclear reactor core. (Only rough estimates are
required.)
2-13. Suppose that the total cross section of rhodium has been measured and the following
values have been obtained for the resonance parameters of a well-isolated resonance at
Eo= 1.26 eV: .a = 5000 b, r = 0.156 eV, and a, = 5.5 b. Plo! the value of the total cross
section for values of the energy between 0.2 and 40 eV. Calculate the thermal
absorption cross section and compare this with the measured value of 156 b. (Assume
that resonance scattering can be neglected.)
2-14. At higher energies, the differential elastic scattering cross section in the CM system
exhibits anisotropy (so-called "p-wave" scattering) of the form
Plot the scattering probability P(Ei-rEr) against final energies Er for this more general
cross section behavior for the three cases a>O, a=O, and a <O. Give a physical
interpretation of your sketches.
2-15. Using the Maxwell-Boltzmann distribution M ( V , T ) , calculate the most probable
energy of the nuclei characterized by such a distribution. Also calculate the average
thermal energy of these nuclei.
2-16. The partial widths of the first resonance in 2 " ~ at 5.49 eV are r,=.029 eV and
r, = ,0018 eV. Plot the Doppler-broadened capture cross section at the temperature of
O0K, 20°C. and 1000°C. [Use the tabulated $(I,%)function.]
2-17. Show that the total area under a Doppler-broadened resonance is essentially inde-
pendent of temperature.
2-18. Using the differential scattering cross section characterizing a proton gas at tempera-
ture T,compute the corresponding macroscopic scattering cross section P,(E). In
particular, determine the behavior of this cross section for low energies E.
2-19. A neutron is absorbed in a "'u nucleus at r=O. Describe a probable life history of the
resulting and its successors on the assumption that it undergoes fission. Give
order of magnitude estimates of characteristic times at which various events occur.
Describe the various particles injected into the system as a result of this fission.
FISSION CHAIN REACTIONS AND NUCLEAR REACYORS--AN INTaODUCIION / 75
such reactions by using the neutron as a chain carrier. It should be apparent that if
we wish to maintain a stable or steady-state chain reaction, that is, one that does
Fission Chain not grow or decay away with time, we must arrange things so that precisely one
Reactions and neutron from each fission will induce another fission event. The remaining fission
neutrons will then either be absorbed in capture reactions or will leak out from the
Nuclear Reactors system. We must design the nuclear reactor to achieve this very delicate balance
between fission reactions and neutron capture and leakage, as we indicate schema-
-an Introduction tically in Figure 3- 1.
We can express this requirement in mathematical form. Since the neutrons play
the central rile in maintaining the fission chain reaction, let us focus our attention
on them for a moment. A given neutron will be "born" in a fission event and will
then usually scatter about the reactor until it meets its eventual "death" in either an
absorption reaction or by leaking out of the reactor. Certain numbers of these
neutrons will be absorbed by fissile or fissionable nuclei and induce further fission,
thereby leading to the birth of new fission neutrons, that is, to a new "generation"
of fission neutrons. Suppose that we could somehow measure the number of
neutrons in two successive fission neutron generations. We would then define the
ratio of these numbers as the multiplication factor k characterizing the chain
reaction
In particular, note that this very simple model of nuclear reactor kinetics agrees with control the reactor power level, for a 0.1% change in the multiplication factor is
our earlier definition of reactor criticality in terms of k (see Figure 3-2). Yet this rather common. Fortunately we have omitted something from this simple model
model also tells us that the growth or decay of the neutron population in a reactor which tends to greatly increase the neutron lifetime 1 and hence T, thereby slowing
obeys an exponential growth law. Such exponential growth is quite commonly down the reactor time response. This is the effect of delayed neutrons on the chain
found in the study of population dynamics. Indeed the study of the "neutron" reaction. However this is a tale for another time, so we will leave our study of
population in a reactor core is mathematically rather similar to the study of reactor kinetics with the promise of returning later to patch up this model in order
biological populations, and hence the terminology of the latter field is frequently to provide a more optimistic picture of nuclear reactor time behavior.
adopted in reactor physics (e.g., generation, birth, life, death, virgin, daughter).
We will later find that the power level of a nuclear reactor is essentially C. A Formal Calculation of k: The Four-Factor Formula
proportional to its neutron population. Hence we can also regard the time behavior
of the reactor power level as being exponential with a time constant or reactor Let us now turn our attention to the calculation of the multiplication factor
period T given by for, say, a pile of uranium that one wishes to make into a nuclear reactor. We
probably should add some coolant to remove fission heat and perhaps some
1 structural material to hold the core together. However we will assume that we can
T==- treat these materials as intimately and homogeneously mixed so that the composi-
k-l '
tion of the reactor is uniform.
In particular it should be noted that as the multiplication factor k approaches Now to calculate k we must determine the possible fate of neutrons in a given
unity, the reactor period T approaches infinity which corresponds to a time- fission generation. Fortunately this is rather easy to do since there are only two
independent neutron population or reactor power level. possible alternative destinies available to the neutron. First it might leak out of the
However suppose that k is not equal to unity. Then how rapidly might we expect reactor and be lost to the chain reaction. If it does not leak out, then it must
the power level of the reactor to change? Suppose, for the sake of illustration, we eventually be absorbed.' This absorption may correspond to a nonproductive
increased k to make the reactor ever so slightly supercritical by an amount of capture event in either the fuel or other materials, or the absorption may induce a
k = 1.001. Since the neutron lifetime in a typical power reactor is about lo-' sec, fission reaction, in which case a new fission neutron generation is produced. We
we find this corresponds to a reactor period of T30.1 sec. Hence in one second the can represent these destinies schematically, as shown below:
power level of the reactor will increase by a factor el0= 22,000. Thus it appears that
the reactor will respond very rapidly to changes in the multiplication factor. In fact,
Leak out of system Absorbed
if a power reactor did indeed respond this rapidly, then it would be difficult to
Fission
neutron
/' or
dAbsorb:l
To make this more formal, suppose we define the probabilities for each of these
possible events as follows:
'Of course, yet a third alternative would be a decay of the neutron into a proton, electron, and
neutrino, but since the half-life for decay of a free neutron is 11.7 minutes, and the typical
FIGURE 3-2. T l w beluvlor ol the number o l neutrons in a reactor. neutron lifetime I in the reactor is less than lo-) sec, we can safely ignore this alternative.
80 / INTRODUCTORY CONCEITS OF NUCLEAR POWER REACI'OR ANALYSIS FISSION CHAIN REACTIONS AND NUCLEAR dEACTORS--AN INTRODUCTION / 81
These latter two conditional probabilities are easily calculated. The conditional where we have recalled that q=v(a;/oh) is the number of fission neutrons
probability for absorption in the fuel PA, can be expressed simply as the ratio of produced per absorption in the fuel. We can now use our definition of the
the macroscopic absorption cross sections for the fuel Xi, and for the fuel plus the multiplication factor k as being the ratio of the number of neutrons in two
rest of the material in the core 2.. (We will usually indicate with a superscript the successive fission generations to write
material to which we are referring. The absence of the superscript will imply that
the macroscopic cross section is the total for all of the materials in the system.)
Thus we can write
The nonleakage probability P,, appearing in Eq. (3-12) and characterizing neutron
leakage from the core is much more difficult to compute. It will require more
elaborate mathematics (and in a realistic calculation, the use of a digital computer),
It should be kept in mind that this expression has been introduced only for the and hence we will defer a discussion of it until later.
situation in which the reactor has a uniform composition. Unfortunately for the As a momentary detour, however, suppose our reactor were of infinite extent.
reactor analyst all modem reactors have nonuniform compositions varying from Then since no neutrons could leak out, we immediately conclude that we must set
point to point (e.g., due to fuel elements, coolant channels, support structure). In the nonleakage probability PNL= I. The corresponding multiplication factor is then
this more general case one can still use Eq. (3-8) if the macroscopic cross sections known as the infinite medium mult$lication factor and denoted by
X, are regarded as spatial averages over the reactor. It should also be noted that we
have not yet specified the neutron energy at which these cross sections are to be
evaluated. Again, we will later find that the cross sections appearing in Eq. (3-8)
must be appropriately averaged over energy, just as they are over space. Now of course no reactor is of infinite size. Nevertheless k , is a useful parameter
It is customary in reactor terminology to refer to this probability as the thermal in reactor analysis since it essentially characterizes the multiplication properties of
utilization of the reactor and denote it by PA,=f. This term arose in the early
the material in the reactor as distinct from the geometry of the reactor core. Of
analysis of thermal reactors in which essentially all fissions in the fuel were induced course since PNL< 1 more generally for a finite reactor from which some neutron
by thermal neutrons. In this case the cross sections in f would be evaluated at leakage can occur, we must have k , > I in order to have any chance of achieving a
thermal neutron energies and would represent the effectiveness of the fuel in critical chain reaction.
competing with other materials in the reactor for the absorption of thermal There are a couple of important modifications that must be introduced into this
neutrons, that is, the effectiveness with which the reactor utilized the thermal simple development in order to understand how the present generation of so-called
neutrons in the fuel. The expression in Eq. (3-8) actually applies to any type of "thermal" reactors works. We must account for the fact that the neutrons in a
nuclear reactor have a distribution of energies. As we saw in Chapter 2, the fission
reactor. However we will fall in line with convention and refer to it as the thermal
neutrons are born at very high energies in the MeV range. However the fission
utilization and denote it by "J."
cross section is largest at very low energies-indeed, at those energies correspond-
The conditional probability for inducing a fission reaction in the fuel can also be
ing to neutrons in thermal equilibrium with the reactor core at a temperature T,
expressed in terms of cross sections. In this case we simply take the ratio of the
e.g., for T=300°C, E = kT=0.05 eV. Hence it is obviously to our advantage to try
fission cross section to that of the absorption cross section (due to both fission and
to slow down, or in the language of reactor physics, "moderate," the fast fission
radiative capture) in the fuel material:
neutrons to take advantage of the fact that slow neutrons are more likely to induce
fission reactions. This can be accomplished rather easily, simply by letting the fast
neutrons collide with light nuclei, thereby losing some of their kinetic energy in
elastic scattering collisions. The lighter the nucleus involved, the more kinetic
energy per collision will be lost on the average by the neutron and hence the more
We are now ready to utilize these probabilities to determine the multiplication effective the slowing down or moderation. In fact the best nucleus to use is
factor k. The general scheme is to play a game of "follow the neutron." Suppose we hydrogen, which is fortunately quite commonly available in the form H,O. Hence
start with N, neutrons present in the reactor in a given fission generation. Then if we just let the fast neutrons rattle around in water for a bit, they will quickly slow
with the help of the above probabilities and our diagram, we can compute the down to the desired thermal energy. In this sense, we refer to water as a neutron
number of neutrons in the next generation as: moderator. Numerous other materials can be used as moderators in nuclear
reactors, and we will discuss these in greater detail later.
The presence of such neutron moderation in a reactor suggests several modifica-
tions to our earlier calculation of the multiplication factor k. Suppose we first
modify our diagram of the various possible neutron destinies to take into account
neutron energy as shown in Figure 3-3. Now since most fissions will be induced by
thermal neutrons, we will regard J and -r) as being evaluated at thermal neutron
82 / INTRODUCTORY CONCEPTS OF NUCLEAR POWER REAtXOR ANALYSIS FISSION CHAIN REACnONS AND NUCLEAR REACTORS--AN INTRODUCnON / 83
The fast fission factor r is usually quite close to unity in a thermal reactor with
typical values ranging between c = 1.03 and t = 1.15.
The second factor we will introduce will characterize the possibility that the
neutron might be absorbed while slowing down from fission to thermal energies.
Since most absorptions occurring during the slowing down process correspond to
resonance capture in heavy nuclei such as "'U, we refer to this factor as the
;----- MeV fission
neutrons
resonance escape probabilify p :
where
-Probability that fast neutron will not leak out
p,=
1 Sub-eV thermal (fast nonleakage)
I neutrons
I lleakage or ---r
L---+- parasitic capture) -Probability that thermal neutron will not leak out
Fission PTNL' (thermal nonleakage).
FIGURE 3-3. Roeeara chwfterlzlq a mutron gemnbn in a lW reactor.
If we now insert these new definitions into our earlier expressions (3-12) and
(3-13), we find that the infinite medium multiplication factor becomes
energies. For example, f would now refer to the ratio of thermal neutron absorp-
tions in the fuel to total thermal neutron absorptions and thereby become more
deserving of its designation as the "thermal" utilization. Similarly, q is now This is known as the four-factor formula. Moreover one now writes
identified as the average number of fission neutrons produced per absorption of a
thermal neutron in the fuel.
Then to account for processes that occur while the neutron is slowing down to which is known, surprisingly enough, as the six-factor formula.
thermal energies, we will introduce two new quantities. We first define a factor that
takes account of the fact that, although most fissions will be induced in fissile
EXAMPLE: To more vividly illustrate these ideas, we can list the values of each
material by thermal neutrons, some fissions will be induced in both fissile and
of the factors in the six-factor formula for a typical thermal reactor: q = 1.65,
fissionable material by fast neutrons. Hence we will scale up our earlier expression
fs0.71, r = 1.02,p=0.87, PFNL=0.97,PT,,=0.99*k,= 1.04*kE 1.00.
for k by a fast fission factor r :
Hence provided we can calculate each of these factors, our criticality condition
Total number of fission neutrons (from both fast and thermal fission) k = l can then be easily checked (in the above example we fudged up the
r z? (3- 14) parameters a bit to yield a critical system). Of course, the calculation of these
Number of fission neutrons from thermal fissions
factors is quite difficult in general. Indeed one cannot really separate the various
84 / INTRODUCTORY CONCEPTS QF NUCLEAR POWER REACXOR ANALYSIS FISSION CHAIN REACTIONS AND NUCLEAR REACTORS-- AN INTRODUCTION / 85
conditional probabilities as was done in these formulas. Instead alternative schemes moderator density, the type of moderator, coolant, and structural materials used, or
based on iterative numerical methods must be used in practice to arrive at a the manner in which reactor multiplication is controlled. One would refer to the
criticality condition. amount of fuel required to achieve a critical chain reaction as the critical mass of
Nevertheless the four-factor and six-factor formulas are quite useful because fuel.
they provide insight into the various mechanisms involved in nuclear fission chain In reality, however, a nuclear reactor is always loaded with much more fuel than
reactions and on rare occasions may actually be of use in making crude estimates is required merely to achieve k = I. For example the LWR is typi~allyloaded with
in nuclear design. They are also useful in illustrating the trends of parameter sufficent fuel to achieve a multiplication of about k = 1.25. This e ~ t r multiplication
a
variation in a core design. is required for several reasons. First if the reactor is to operate at power for a
For example, although and c are essentially fixed once the fuel has been period of time, one must provide enough excess fuel to compensate for those fuel
chosen, the thermal utilization f and resonance escape probability p can be varied nuclei destroyed in fission reactions during the power production. Since most
considerably by changing the ratio of fuel density to moderator density. All of contemporary reactors are run roughly one year between refueling, a sizable
these parameters can be varied by using a heterogeneous lattice of fuel elements amount of excess fuel is needed to compensate for fuel burnup. A second
surrounded by moderator rather than a uniform, homogeneous mixture of fuel and motivation arises from the fact that the multiplication of a reactor tends to
moderator. decrease as the reactor power level and temperature increase from ambient levels to
One can also vary the nonleakage probabilities by simply making the reactor operating levels. Additional multiplication is needed to compensate for this effect.
core larger, or surrounding the reactor by a material with large scattering cross Finally one must include enough extra multiplication to allow for reactor power
section so that some of the neutrons leaking out will be scattered back into the level changes. For example, we have seen that if we wish to increase the reactor
reactor. Actually when leakage is changed, there will be some change in the power level, we must temporarily adjust k to a value slightly greater than 1 so that
parameters in the four-factor formula as well since these are actually averages over the reactor is supercritical. The reactor can then be returned to critical when the
the various neutron energies in the reactor, and this distribution of enerses will desired power level has been reached.
vary with the amount of leakage. Such considerations have given rise to a Of course when this excess multiplication is not being used, some mechanism has
somewhat different notation for the multiplication factor characterizing a finite to be provided to cancel it out to achieve reactor criticality. This is the function of
system which is occasionally referred to as the effecriw multiplicarion factor and reactor control mechanisms. Such control is usually achieved by introducing into
denoted by kd, the reactor core materials characterized by large absorption cross sections. They
will then tend to eat up the excess neutrons produced in the chain reaction. In
terms of our six-factor formula such absorbing materials lower the value of the
thermal utilization f, since they compete with the fuel for neutron absorption. A
There are other prescriptions for defining the multiplication factor. In particular we variety of types of reactor control are used in power reactors. For example, the
will introduce one of these schemes later when we consider the analytical treatment neutron absorber might be fabricated into rods which can then be inserted into or
of the neutron energy dependence in more detail. However for now we will withdrawn from the reactor at will to vary multiplication. Sometimes the absorber
continue to regard the multiplication factor as the ratio between either the number is fabricated directly into the fuel itself. Or it may be dissolved in the reactor
of neutrons in two successive fission generations or the neutron production and coolant. When such control absorbers are used to hold down the excess multiplica-
loss rates in the reactor. tion introduced to compensate for fuel burnup, one refers to them as shim control.
We can use the six-factor formula for the multiplication factor to gain a bit more They may also be used to force the reactor subcritical in the case of an emergency;
insight into the goals of reactor design and operation. There are several ways to then they are known as scram control. Finally they may just be used to regulate the
adjust k in the initial design of the reactor. One could first regard the size of the power level of the reactor; then they are referred to as maneuwring control
reactor as the design variable. Since the ratio of surface area to volume decreases as elements.
the reactor geometry is enlarged, one can control the relative importance of the The ease with which such control elements can control the fission chain reaction
leakage factors by adjusting the reactor size., For a given wre composition (with k, will depend on how rapidly the reactor responds to variations in multiplication.
greater than 1, of course) there will be a certain critical size at which k = I. An Since fuel burnup occurs over very long periods of time (typically weeks or
alternative way to achieve the same reduction in leakage is to surround the reactor months), a rapid response of shim control is not required-which is fortunate,
with a scattering material that acts as a neutron reflector. Most thermal reactor because rather large amounts of multiplication must be manipulated (typically
cores are so large that leakage represents a rather small loss mechanism (typically changes of 10-2096 in k). The normal power variations in the reactor are due to
about 3% of the neutrons leak out from the core in large thermal reactors). much smaller changes in multiplication (<0.1%) and are characterized by essen-
Usually the core size and geometry for a power reactor are dictated by thermal tially the reactor period T which in turn is proportional to the neutron lifetime I.
considerations, for instance, the size of the core necessary to produce a given power However we saw earlier that the lifetime of prompt fission neutrons was quite
output while being provided with sufficient cooling so that the temperature of the short, typically about sec. The effective neutron lifetime is greatly increased
reactor materials will not become excessively high. The primary design variable at by the presence of delayed neutrons, however. We recall that about 0.7% of the
the disposal of the nuclear engineer is the core composition. In particular he can neutrons produced in fission are delayed anywhere from 0.6 to 80 sec since they
vary the composition (enrichment) and shape of the fuel, the ratio of fuel to arise from fission product radioactive decay. Hence the effective neutron lifetime is
86 / INTRODUCTORY CONCEPTS OF NUCLEAR WWER REACTOR ANALYSIS FISSION CHAIN REACITONS AND NUCLEAR REACTORS-AN INTRODUCllON / 87
actually the average of the prompt neutron lifetime and the average decay time of will leak out or be absorbed in parasitic capture) while one neutron will be needed
these delayed neutrons, properly weighted, of course, by their relative yield frac- to replace the consumed fissile nucleus by converting a fertile into a fissile nucleus.
tions. When this is taken into account, one finds that the effective neutron lifetime If we return to Figure 2-25 we can see that the only attractive breeding cycle for
is almost two orders of magnitude longer, Id,-- lo-' seconds. Hence a multiplica- low-energy (i.e., thermal) neutrons would involve 2 3 3 ~that , is, the 2 3 2 ~ h / 2 3 3 ~
tion of 0.1% would now correspond to a reactor period of T = 10 seconds, well process. To breed using U 8 U / 2 3 9requires
~~ that we use fast neutrons with energies
within the control capability of a reactor control system. greater than 100 keV. And of course this is the motivation behind the development
of the fast breeder reactor.
D. Conversion and Breeding At this point, it is useful to digress a bit and discuss the average energy of the
neutrons sustaining the chain reaction in various types of nuclear reactors. As we
If we recall our earlier expression for the multiplication factor k in Eq. (3-12), have seen, the energies of neutrons in a reactor span an enormous range, from 10
it is evident that since the thermal utilization f and the nonleakage probability P,, MeV (usually the maximum energy of fission neutrons) down to as low as eV
are both less than I, we require TJ to be substantially greater than 1 if a critical after having suffered a number of scattering collisions with nuclei and slowing
fission chain reaction is to be possible. Fortunately as we can see from Figure 2-25, down. Furthermore the neutron cross sections depend sensitively on the neutron
this condition is not only satisfied, but in fact for many energies one finds that energy. As the examples in Chapter 2 indicated, the general trend is for cross
v > 2. Hence we in fact appear to have an extra neutron. This "bonus" neutron can sections to decrease with increasing energies. This feature is particularly true of
be put to good use if we recall that certain fertile isotopes can be transmuted into absorption cross sections such as capture or fission.
fissile material via neutron capture. In particular, UBU can be transmuted into The fact that the fission cross section of is largest at low energies implies that it is
2 3 9 while
~ ~ 13'Th can be transmuted into 233U. Hence if we load the core of a easiest to maintain a fission chain reaction using slow neutrons. Hence early
reactor with such fertile material, we can use the extra neutron to produce a new nuclear reactors used low mass number materials such as water or graphite to slow
fissile fuel material. This process is frequently referred to as conwrsion, and nuclear down or moderate the fast fission neutrons. Such moderating materials slow the
reactors whose principal job is to produce 2 3 9 Por ~ U3U are known as conwrter neutrons down to energies comparable to the thermal energies of the nuclei in the
reactors. reactor core. Reactors characterized by an average neutron energy comparable to
Actually all modem power reactors are converter reactors in a sense, although such thermal energies are referred to as thermal reactors. Such reactors require the
this is not their primary function, since they contain substantial amounts of 238U minimum amount of fissile material for fueling and are the simplest reactor types
which will be transmuted into 2 3 9 Pvia
~ neutron capture during normal operation. to build and operate. Most nuclear power plants in this country and abroad utilize
For example, a LWR will contain a fuel mixture of roughly 3% U5U and 97% U8U thermal reactors.
in a freshly loaded core. After a standard operating cycle (three years), this However we have also seen that there is a very definite advantage in keeping the
fuel will contain roughly 1% 235Uand I% 2 3 9 Pwhich
~ can then be separated out of neutron energy high, since the number of neutrons emitted per neutron absorbed in
the spent fuel and refabricated into fresh fuel elements for reloading (so-called the fuel TJ is largest for fast neutrons. Hence one can use the "extra" neutrons
"plutonium recycling"). available in a fission chain reaction maintained by fast neutrons to convert or
These considerations suggest that it might in fact be possible to fuel a reactor breed new fuel. However since a, is smaller, one also needs much more fuel to
with 'I9pu and 13'U and then produce directly the fuel (U9Pu) needed for future sustain the chain reaction. Furthermore to keep the neutron energy high, one wants
operation. Indeed it might even be possible to produce more 2 3 9 P ~than is to utilize only high mass-number materials in the core to keep neutron slowing
burned-that is, to "breed" new fuel. This is the essential idea behind the concept down to a minimum. Such reactors characterized by average neutron energies
of a breeder reactor. above 100 keV are known as fast reactors. It is felt by many that fast reactors will
To discuss this concept in more detail, it is useful to define the conwrsion ratio eventually replace the current generation of thermal power reactors because of
their ability to breed fuel.
Average rate of fissile atom production To make some of these ideas a bit more precise, we have compared the
CRs Average rate of fissile atom consumption important nuclear parameters v, TJ,and of for typical nuclear fuels at energies
characterizing both thermal and fast reactors. (To be more precise, these quantities
This quantity is also referred to as the breeding ratio (BR) if it is greater than one. have been calculated by averaging the energy-dependent nuclear parameters v(E),
If we have conversion then, consuming N atoms of fuel during reactor operation q(E), and a,(E) over the neutron energy distributions found in typical LWRs and
will yield CR.N atoms of the new fissile isotopes. For example, most modern LMFBRs.) One should first note that the fission cross sections in fast reactors are
LWRs are characterized by a conversion ratio of CRs0.6. By way of contrast, some two orders of magnitude lower than those in thermal reactors. Hence even
HTGRs are characterized by somewhat higher conversion ratios CR20.8 and though fast reactors exhibit considerably higher conversion ratios (typically, CR
hence are sometimes referred to as advanced converter reactors. = BR- 1.2-1.5) due to a larger value of TJ,their fissile inventory requirements may
For breeding to occur we require that the conversion ratio be greater than unity, run as much as several times those required by thermal reactors just to maintain a
CR = BR > 1. Of course for this to happen we must have q > 2 since slightly more critical chain reaction with fast neutrons. This table also indicates that while 235U
than one fission neutron is needed to maintain the chain reaction (some neutrons will yield a slightly higher conversion ratio than 'j9pu in thermal reactors, the use
88 / INTRODUCTORY CONCEPTS OF NUCLEAR POWER REACTOR ANALYSIS
of 2 3 9 does
~ ~ exhibit a sizable advantage in fast reactors since the capture-to-fission
ratio ~y~~ falls off quite markedly for large neutron energies. In Table 3-1 we have
only indicated the nuclear fission properties of 238Uin fast reactors, since this
isotope is fissionable and hence contributes only a modest fraction of the fissions
occurring in thermal reactors (-2-5%) in contrast to its rather large contribution in
fast reactors (-20%).
temperature cooling water which is usually obtained from artificial cooling ponds or
cooling towers.
This is of course a very oversimplified description of the major components of a
power plant, but it does serve to illustrate that these components are quite similar
for both nuclear and fossil-fueled stations. Actually as far as the steam cycle itself
is concerned, the primary difference between the two types of plant is that the
fossil-fueled boiler supplies slightly higher temperature, higher pressure steam,
thereby reducing the design requirements on the turbine (although it should be
mentioned that more advanced reactor types such as the HTGR supply system
steam at conditions quite comparable to those of modern fossil-fueled units). Of
course there are numerous other differences in the various subsystems of the plants,
as well as in their operation. However the major features of the plants, aside from
their steam supply system, are quite similar.
Steam - il fi
LMFBR
coolant loop must be utilized to isolate the steam generator from the very high
induced radioactivity of the primary coolant loop passing through the reactor.
The most common coolant used in power reactors today is ordinary water, which
serves as both coolant and moderating material in the reactor. There are two major
types of LWR: pressurized water reactors (PWR) and boiling water reactors
(BWR). In a PWR the primary coolant is water maintained under very high
I pressure (-155 bar) to allow high coolant temperatures without steam formation
m I
Multiple-stags turbine
within the reactor. The heat transported out of the reactor core by the primary
I- 1
coolant is then transferred to a secondary loop containing the "working fluid" by a
Condenser steam generator. Such systems typically contain from two to four primary coolant
loops and associated steam generators.
In a BWR, the primary coolant water is maintained at a sufficiently low pressure
Feedwater pump
(-70 bar) for appreciable boiling and steam formation to occur within the reactor
Transformer
core itself. In this sense the reactor itself serves as the steam generator, thereby
Faadwater eliminating the need for a secondary loop and heat exchanger. In both the PWR
Feedwater heaters
Balance of plant (1 ([
\7 and BWR, the nuclear reactor itself and the primary coolant are contained in a
large steel pressure wssel designed to accomodate the high coolant pressures and
temperatures. In a PWR, this pressure vessel must be fabricated with thick steel
walls to contain the very high primary coolant pressures. By way of contrast, the
FIGURE 3-5. Comparison of different steam-supply systems. BWR pressure vessel need not be so thick, but must be much larger to contain both
92 / INTRODUCTORY CONCEPE3 OF NUCLEAR W W E R REACTOR ANALYSIS FISSION CHAIN REA(TI?ONS AND NUCLEAR REACTORS-AN INTRODU(XI0N / 93
(4) Moderator: Material of low mass number which is inserted into the
reactor to slow down or moderate neutrons via scattering collisions.
End cap, Typical moderators include light water, heavy water, graphite, and
beryllium.
(5) Coolant: A fluid which circulates through the reactor removing fission
heat. The coolant can be either liquid, such as water or sodium, or
gaseous, such as helium or carbon dioxide. It may also serve a dual role
as both coolant and moderator, such as in the LWR.
(6) Coolant channel: One of the many channels through which coolant flows
in the fuel lattice. This may be an actual cylindrical channel in the fuel
assembly, as in the HTGR, or an equivalent channel associated with a
single fuel rod, as in a LWR.
(7) Structure: The geometry and integrity of the reactor core is maintained
by structural elements such as support plates, spacer grids, or the
metallic tubes used to clad the fuel in some reactor designs. The
structural materials may also serve a dual role by moderating neutrons
Pellet such as the graphite in an HTGR.
( 8 ) Control elements: Absorbing material inserted into the reactor to control
core multiplication. Although most commonly regarded as movable rods
of absorber, control elements may also consist of fixed absorbers or
absorbing materials dissolved in the coolant. Common absorbing
materials include boron, cadmium, gadolinium, and hafnium.
(9) Reactor core: The total array of fuel, moderator, and control elements.
(10) Reactor blanket: In a breeder or high conversion reactor the core is
usually surrounded by a blanket of fertile material that more effectively
utilizes the neutrons leaking out of the core.
(11) Reflector: A material characterized by a low absorption cross section
Fuel ekment
" -
Fuel assembly FIGURE 3-7. Fwd ekaent; Inel ssembly.
used to surround the core in order to reflect or scatter leaking neutrons
back into the core.
(12) Shielding: The reactor is an intense source of radiation. Not only must
(I) Fuel: Any fissionable material. This can be either fissile material such as operating personnel and the public be shielded from this radiation, but
"%, 235U,139P~, or 2 4 1 Por
~ fissionable material such as 232Th,2 3 8 ~or, reactor components must as well be protected. Hence absorbing material
is introduced to attenuate both neutron and gamma radiation. Thermal
2 4 0 ~ Most
~ . modern power reactors utilize this fuel in a ceramic form-
shielding is used to attenuate the emergent core radiation to levels that
either as an oxide such as UO,, a carbide such as UC, or a nitride, UN.
do not result in significant heat generation and hence damage in reactor
(2) Fuel element: The smallest sealed unit of fuel. In an LWR or LMFBR
the fuel element is a metal tube containing ceramic pellets of fuel (such components. Biological shielding reduces the radiation still further to
acceptable levels for operating personnel.
as U02). (See Figure 3-7.) In an HTGR the fuel element can be regarded
(13) Support structure: The support plates that serve to maintain the core
as either a tiny (300 pm diameter) particle of uranium carbide coated
with pyrolytic graphite layers, or as a cylindrical fuel pin composed of geometry.
these fuel particles bound together with a graphite binder. (14) Reactor pressure vessel: The high pressure containment for reactor and
(3) Fuel assembly or bundle: The smallest unit combining fuel elements into associated primary coolant system.
an assembly. For example, in a LWR the fuel assembly is composed of It is also useful to introduce at this point several quantities which are used to
several hundred fuel elements fastened together at top and bottom with describe reactor performance. The units in which these quantities are usually
coolant nozzle plates and with several spring clip assemblies along the expressed are denoted in brackets.
length of the fuel (see Figure 3-7). In an HTGR the fuel assembly is a
hexagonal block of graphite with holes into which the cylindrical fuel (I) Reactor thermal power [MWt]: The total heat produced in the reactor
pins are inserted. Fuel is usually loaded into a reactor core or replaced core.
one fuel assembly at a time. A typical power reactor core will contain (2) Plant electrical output [MWe]: Net electrical power generated by the
hundreds of such fuel assemblies. plant.
% / INTRODUCTORY CONCEPTS OF NUCLEAR POWER REACTOR ANALYSIS FISSION CHAIN REACTIONS AND NUCLEAR REACTORS-AN INTRODUCnON / 97
Plant electrical output straints imposed on the reactor operation. The nuclear analysis and design of a
(3) Net plant efficiency [%I:
Reactor thermal power ' reactor core is highly dependent on other areas of core design, including thermal-
Total energy generated over time period hydraulic design, structural analysis, economic performance, and so on. The
(4) Plant capacity factor [%I: criteria for a design effort are quite varied, encompassing considerations of per-
(Plant rating) X (time)
formance, reliability, economics, and safety. These criteria are frequently contra-
Average plant electrical power level dictory in nature, and hence require optimization.
( 5 ) Plant load factor [%I:
Peak power level The complete nuclear design of a given core configuration is performed many
times, initially to survey design parameters, identify constraints, then to refme the
Integrated electrical energy output capacity
(6) Plant aoailabili~yfactor [%I: design while interacting with other facets of the plant design, and finally, to
Total rated energy capacity for period ' establish a reference design that provides a calculational base against which
optimization calculations can be compared.'
Reactor thermal power This design process is very similar to that utilized in other fields of engineering.
(7) Core power density [kW/liter]: One first must attempt to define the various design constraints that include
Total core volume '
considerations of system performance in terms of both system reliability and
( 8 ) Lineorpower density [kW/m]: Thermal heat generated per unit length of economic performance and safety criteria. Next a preliminary design is proposed,
coolant channel. drawing on available information such as plants already in operation, experimental
Reactor thermal power mockups, and frequently, old-fashioned intuition. Such a design includes a set of
(9) Specific power [kW/kg]: - Total mass of fissionable material ' specifications involving quantities such as fuel enrichment, coolant flow rates and
temperatures, core configurations, reload patterns, and so on. A detailed analysis
(10) Fissile loading [kg]: Total mass of fissionable material. of this preliminary design is then performed in order to evaluate its predicted
Mass of fissile material performance and ascertain whether it conforms to the constraints imposed on the
(I I) Fuel enrichment [%I:
Mass of fissile and fertile material ' system. For example, one would want to calculate the core power and temperature
(12) Fuel burnup [Megawatt-days/metric ton uranium = MWD/TU]: distribution, the pressure drop of the coolant as it passes through the core,
coolant-flow conditions, and the fuel lifetime. When possible, these calculations are
Energy generated in fuel during core residence compared against experiments in order to validate the computational models used.
Total mass of fuel A detailed evaluation of the preliminary design will then lead to more detailed
Fuel burnup designs and analyses as one attempts to optimize the tradeoff between system
(13) Fuel residence time: performance and design constraints. As a final design is approached, one attempts
(Specific power) X (capacity factor) '
to define detailed system specifications.
These are the more common terms used in characterizing nuclear plant perfor- The above procedures emphasize the importance of adequate models of a
mance. We will introduce other more specific concepts and terminology later as we nuclear reactor in order to carry out the required parameter and optimization
develop the more detailed theory of nuclear reactor behavior. studies. These models must be realistic since nuclear reactors are fa1 too expensive
to be built without detailed and accurate design information. Unfortunately any
calculation sufficiently realistic to be of use in reactor design is far too complex to
111. NUCLEAR REACTOR DESIGN be carried out by hand. Hence the digital computer plays a very key role in nuclear
reactor d e ~ i g n . ~
A. General Design Functions of the Nuclear Engineer A key task of the nuclear reactor engineer is to develop models of nuclear
reactors that can then be analyzed on the computer. Such models result in large
The design of a large nuclear power plant is an enormously complex task and
involves the coordination of a remarkably diverse range of disciplines. Each major computer programs or "codes" which can then be used by other nuclear engineers
component of the plant requires a separate and distinct design analysis and is in reactor design. Most of our emphasis in this text is on learning how to synthesize
usually the responsibility of a specific engineering design team. For example, the such approximate models of reactor behavior and then cast them in a form suitable
design of the reactor pressure vessel or steam generators is usually performed by for reactor design. In the language of nuclear engineering, then, this text should be
the reactor supplier, while the turbogenerator and switchgear design is the re- regarded as a primer on nuclear methods deoelopment.
sponsibility of the electrical equipment manufacturer. The coordination among A word of caution should be inserted here, however. In the early days of reactor
these different design projects is extremely important, however, since the designs development it was hoped that one would eventually be able to accurately model
frequently interact to a very high degree. nuclear reactor behavior utilizing only fundamental principles and measured
The primary responsibility for the nuclear design of the reactor core rests with nuclear data. However over the past three decades of reactor development ex-
the nuclear engineer. This design must be accomplished within numerous con- perience it has become apparent that the accuracy of nuclear data and computa-
tionally feasible analytical methods are simply not sufficient to allow this6 Instead
98 / INTRODUCTORY CONCEPTS OF NUCLEAR POWER REACTOR ANALYSIS FISSION CHAIN REACTIONS AND NUCLEAR REACTORSAN INTRODUCTION / 99
nuclear reactor analysis has relied heavily on that very basic ingredient utilized in Our theoretical approach is to begin with an essentially exact description of the
most other areas of engineering design known as the "enlightened fudge." That is, neutron density in the reactor based on the so-called neutron transport equation.
most nuclear analysis methods or computer codes contain empirical parameters This equation, while relatively easy to derive, is extremely difficult to solve, and
that have been adjusted or calibrated by comparing the predictions of the methods hence we will be concerned with developing various approximations to it that !end
with actual experimental measurements. While such empirical input is usually very
themselves more readily to practical application. We begin our actual study of
successful in yielding accurate nuclear design information with a minimum amount
of effort, the novice nuclear engineer should approach any existing nuclear analyti- nuclear reactor theory by using the simplest such approximation, that in which the
cal method or computer code with a high degree of skepticism, since methods neutron energy dependence is neglected by assuming all neutrons to be
calibrated to work well for one range of parameters may fail miserably when characterized as having a single speed, and describing their transport from point to
applied to new situations in which only limited experience is available. point as a simple diffusion process. This very simple model suffices to develop most
The intimate relation between computers and reactor design cannot be over- of the concepts, as well as to illustrate most of the practical computational
stressed. It is almost impossible for the present-day nuclear engineer to function techniques used in more detailed reactor analysis.
without a reasonable background in computer techniques (both in programing and We next develop a more sophisticated model of the neutron density behavior
numerical analysis). Nevertheless the increasingly heavy reliance of the nuclear based on breaking up the range of neutron energies into intervals or "groups" and
reactor industry on elaborate computational models of reactor performance makes then describing the diffusion of neutrons in each of these groups separately,
it even more imperative that the nuclear engineer possess a very thorough accounting for the transfer of neutrons between groups caused by scattering. Such
background in the fundamental physical and mathematical concepts underlying muftigroup diffusion models are the principal tools used in modern reactor analysis,
these models, as well as a healthy dose of skepticism when he attempts to utilke and we consider them in some detail.
their predictions in reactor analysis. In the final section of the book we illustrate these models by applying them to
analyze several typical problems encountered in nuclear reactor design. In particu-
lar, we explore the relation between such nuclear analysis methods and the other
B. Some Concluding Remarks types of analysis required in nuclear reactor core design.
In these last three chapters we have attempted to introduce several simple but
important concepts involved in nuclear fission chain reactions. We have also
provided a brief overview of nuclear reactor systems and the function of the
nuclear engineer in the design of such systems. With this background we now turn REFERENCES
our attention to a development of the theory underlying the nuclear analysis of
fission reactors. We have seen that the neutron plays the central role as the chain 1. A. M. Weinberg and E. P. Wigner, The Physical Theory of Neutron Chain Reactors, The
carrier perpetuating the chain reaction. The key problem of reactor theory, then, is University of Chicago Press (1958).
to determine the distribution of neutrons in a reactor core. This will not only allow 2. S. Glasstone and A. Sesonske, Nuclear Reactor Engineering, 2nd Ed., Van Nostrand,
one to study the chain reaction process itself, but, as we will later find, since the Princeton, N.J. (1975).
neutron density is proportional to the rate at which fission reactions occur and 3. M. M. El-Wakil, Nuclear Energy Conwrsion, Intext, Scranton (1971).
hence proportional to the core power density, the neutron density is also the key to 4. There are a number of other useful sources on nuclear power systems. Each of the major
the subsequent thermal and mechanical analysis of the reactor. As we have suppliers of NSSS prepare detailed systems descriptions. For example, a very informative
mentioned earlier, there are essentially two aspects to this problem. reference is
One must study the interaction of neutrons with matter-specifically, with the Systems Summary of a Westinghouse Pressurized Water Reactor Nuclear Power Plant,
nuclei that make up the matter. This amounts to either experimental or theoretical Westinghouse Electric Corporation (1971).
The most detailed descriptions of modern nuclear power plants can be found in the
determination of the probabilities that various neutron-nuclear interactions will multivolume set of Preliminary Safety Analysis Reports (PSARs) or Final Safety Analysis
occur-that is, a determination of the appropriate neutron-nuclear cross sections. Reports (FSARs) prepared for each nuclear plant. Of particular interest are the standard
This, however, is not the principal concern of the theory we will develop in this text safety analysis reports prepared for each of the major NSSS types. For example:
but is more properly the domain of the nuclear physicist. Hence we tend to take Babcock and Wilcox Standard Nuclear Steam System, B-SAR-241 (1974).
microsco?ic cross section data as given (in a form to be discussed later), and turn BWR/6 Standard Safety Analysis Report, General Electric Company (1973).
our attention instead to the manner in which these data are utilized in nuclear CESSAR, Combustion Engineering Standard Safety Analysis Report, System 80,
reactor analysis. Combustion Engineering (1973).
Of comparable importance is the study of the transport or diffusion of neutrons GASSAR 6, General Atomic Standard Safety Analysis Report, GA-A13200 (1975).
within a nuclear reactor core as they stream around inside the core, suffering Clinch River Breeder Reactor Plant, Reference Design Report, Westinghouse Electric
collisions with nuclei, occasionally being absorbed, inducing fission reactions, or Corporation, 1974; PSAR (1975).
leaking out through the surface of the core. It is this latter study that will allow us RESAR-3, Reference Safety Analysis Report, Westinghouse Nuclear Energy Systems
to develop models for calculating the distribution of neutrons within the reactor (1973).
5. A. Sesonske, Nuclear Power Plant Design Analysis, USAEC Document TID-26241 (1973).
core. 6. J. Chernick, Reactor Technol. 13. 368 (1971).
100 / INTRODUCTORY CONCEPTS OF NUCLEAR POWER REACM)R ANALYSIS
PROBLEMS
3-1 What is the maximum value of the multiplication factor that can be achieved in any
conceivable reactor design?
3-2 Using the alternative definition of the multiplication factor based on the concept of
neutron balance, repeat the derivation of the six-factor formula.
3-3 A spherical reactor wmposed of "'U metal is operating in a critical steady state.
The One-Speed
Discuss what probably happens to the multiplication of the reactor and why, if the
system is modified in the following ways (treat each modification separately, not
cumulatively): (a) the reactor is rapidly compressed to one-half its original volume, (b)
Diffusion Model
a large, fat reactor operator accidentally sits on the reactor, squashing it into an
ellipsoidal shape, (c) a thick sheet of cadmium is wrapped around the outside of the
reactor, (d) the reactor is suddenly immersed in a large container of water, (e) a source
of a
of neutrons is placed near the reactor, (I) another identical reactor is placed a short
distance from the original reactor, and (g) one simply leaves the reactor alone for a
period of time.
Nuclear Reactor
3 4 One defines the doubling time for a breeder reactor as the amount of time required for
the original fissile loading of the reactor to double. Find an expression for the doubling
time td in terms of: (a) the original fissile loading M,, (b) the power level of the reactor
P a wfF,where Ff is the fission rate occurring in the reactor core, (c) and the breeding
ratio BR.
3-5 A detailed comparison of typical power reactor core parameters is given in Appendix
H. Choose one of the reactor types in this Appendix and perform the following
calculations: (a) verify that the average linear power density, power density, and
specific power given in the table are consistent with the wre volume, thermal power
rating, and fuel loading, @) determine the discharge fuel burnup when the capacity
factor of the nuclear unit is 80% and the fuel residence time is three years, (c)
determine a range for core height and wre diameter and sketch a wre cross-section for
one array of assemblies using the tabulated core data.
3-6 Calculate and plot k, as a function of enrichment from 0.7% "5U to 1005%U5U.Use
the thermal cross section data of Appendix A and assume p = c = I.
3-7 Derive a relationship between the waste heat rejected from a plant of a given output
and the thermal efficiency of the plant. (Several years ago such waste heat was referred
to as "thermal pollution." In a countermove, several of the more optimistic spokesmen
for the nuclear power industry coined the phrase "thermal enrichment.") Using this
expression, estimate the waste heat rejected by: (a) a modern fossil-fuel plant, (b) a
LWR plant, (c) a HTGR plant, (d) an LMFBR plant, and (e) a fusion reactor plant,
assuming that all of these plants are rated at 1000 MWe. Treat the efficiency of the
plant as that for an ideal (Carnot) heat engine.
3-8 Consider an infinitely large homogeneous mixture of 23'U and a moderating material.
Determine the ratio of fuel-to-moderator density that will render this system critical for
the following moderators: (a) graphite, (b) beryllium, (c) water (H,O), and (d) heavy
water ( 4 0 ) . Use the thermal cross section data given in Appendix A.
3-9 Modify the simple description of the time behavior of the neutron population in a
reactor given by Eq. (3-5) to account for the presence of a source in the reactor
producing So neutrons per second. In particular, determine the time behavior of the
neutron population for each of the three cases: k < 1, k = I, and k > 1.
Neutron Transport
We now turn our attention to the central problem of nuclear reactor theory, the
determination of the distribution of neutrons in the reactor. For it is the neutron
distribution that determines the rate at which various nuclear reactions occur
within the reactor. Furthermore by studying the behavior of the neutron population
we will be able to infer the stability of the fission chain reaction. To determine the
distribution of neutrons in the reactor we must investigate the process of neutron
transport, that is, the motion of the neutrons as they stream about the reactor core,
frequently scattering off of atomic nuclei and eventually either being absorbed or
leaking out of the reactor. Most reactor studies treat the neutron motion as a
diffusion process. In effect one assumes that neutrons tend to diffuse from regions
of high neutron density to low neutron density, much as heat diffuses from regions
of high to low temperature, or even more analogously, as one gas of molecules
(corresponding to the neutrons) would diffuse through another (the nuclei) to
reduce spatial variations in concentration.
Unfortunately, however, while the treatment of thermal conduction and gaseous
diffusion as diffusion processes is usually found to be quite accurate, the treat-
ment of neutron transport as a diffusion process has only limited validity. The
reason for this failure is easily understood when it is noted that in most diffusion
processes the diffusing particles are characterized by very frequent collisions that
give rise to very irregular, almost random, zigzag trajectories. However, we have
seen that the cross section for neutron-nuclear collisions is quite small (about lo-"
cm2). Hence neutrons tend to stream relatively large distances between interactions
(recall that the mean free path characterizing fast neutrons is typically on the order
of centimeters). Furthermore, the dimensions characterizing changes in reactor core
composition are usually comparable to a neutron mfp (e.g., a reactor fuel pin is
typically about 1 cm in diameter).
103
161 / THE ONESPEED DIFFUSION MODEL OF A NUCLEAR REACTOR NEUTRON TRANSPORT / 105
Hence one frequently requires a more accurate description of neutron transport chapter since we recognize that many nonnuclear engineers reading this text may
that takes into account the relatively long neutron mfp and neutron streaming. not really need (or care) to understand the limitations and ranges of validity of
Such a description has been borrowed from the kinetic theory of rarefied gases nuclear analysis methods. Therefore if the reader is faint of heart in confronting the
(which are also characterized by long mfp)--more precisely, the kinetic theory of transport equation, and strong in faith in being able to accept the rather heuristic
gas mixtures. The fundamental equation describing dilute gases was first proposed arguments necessary to develop these approximate theories (e.g., the neutron
more than one century ago by Boltzmann, and even today the Boltzmann equation diffusion equation) without recourse to the transport equation, or perhaps just
remains the principal tool of the gas dynamicist.' Its counterpart for the neutron disinterested, he can proceed immediately to the development of neutron diffusion
"gas," the so-called neutron transport equation, is far younger (less than 40 years theory in the next chapter.
old), far simpler (e.g., it is a linear equation in contrast to the Boltzmann equation, The more formal discussion in this chapter will also serve to introduce the
which is nonlinear), but usually strikes far more terror in the hearts of fledgling standard numerical approximation schemes used to analyze the neutron behavior
nuclear engineers who are intimidated by its frightening reputation within the in nuclear reactors. As in other areas of physical analysis, we will find that the
nuclear reactor community. Neutron transport theory has come to be associated usual maxim applies: that the "brute force" numerical approach (i.e., discretize
with a hideous plethora of impenetrable mathematics, unwieldy formulas, and everything in sight and slap in on a computer) is conceptually the simplest
(eventually) the expenditure of enormous amounts of money on computer number- approach to understand and wmputationally the most expensiw calculation to
crunching. perform. The more elegant approximate methods require far less computational
This is most unfortunate because the neutron transport equation is much simpler effort but far more in the way of mental gymnastics in order to understand the
to derive (requiring only the concept of neutron conservation plus a bit of significance and reliability of their predictions.
vector calculus) and to understand than the neutron diffusion equation that we
shall utilize in most of our development of reactor analysis. It is also a far more
fundamental and exact description of the neutron population in a reactor-indeed, I. INTRODUCTORY CONCEPTS
it is the fundamental cornerstone on which all of the various approximate methods
used in nuclear reactor analysis are based. A. Neutron Density and Flux
It does have one major drawback, however. It is usually very difficult to solve
the transport equation for any but the simplest modeled problems (and even these Our ultimate goal is to determine the distribution of neutrons in a nuclear
reactor core. This requires accounting for the neutron motion about the wre and
require an inordinate amount of analytical work). However that is quite all right, neutron interactions with nuclei in the wre. We will begin by defining the neutron
since it is not our intent to attack the transport equation head on. Rather the job of density N (r, t) at any point r in the reactor core by
the reactor analyst is to develop suitable (i.e., calculationally feasible and accurate)
approximations to it. Usually, however, only by comparing these various approxi- N (r, t) d3r=expected number of neutrons in
mate theories to the transport equation from which they originated can one really d3r about r at a time t. (4- 1)
assess their range of validity.
There is another reason for including an introduction to the neutron transport
equation in even an elementary discussion of nuclear reactor analysis. Although The word "expected" has been inserted into this definition to indicate that this will
neutron diffusion theory is usually found adequate for reactor applications, it owes be a statistical theory in which only mean or average values are calculated. (The
its accuracy to various schemes that have been developed to "patch it up" using actual neutron density one would obtain from a series of measurements would
results from more accurate transport equation solutions. For example, we will find fluctuate about this mean value, of course.) The neutron density N (r,f) is of
that the neutron diffusion equation is quite invalid near the boundary of a reactor, interest because it allows us to calculate the rate at which nuclear reactions are
or near a highly absorbing material such as a fuel rod or a control element. occurring at any point in the reactor. To understand this, let us suppose for
Nevertheless we can continue to use diffusion theory to describe the reactor convenience that all the neutrons in the reactor have the same speed u. Now recall
provided we fudge it a bit by inserting so-called "transport corrections" into the that one can express the frequency with which a neutron will experience a given
boundary conditions accompanying the diffusion equation. neutron-nuclear reaction in terms of the macroscopic cross section characterizing
So hopefully we have made a case for our inclusion of a very introductory that reaction Z and the neutron speed u as
discussion of neutron transport theory within an elementary text. We would
caution the reader not to be intimidated by the notation or the apparent oX = interaction frequency. (4-2)
strangeness of the equation we will develop. He should find it rather easy to
understand the derivation and interpretation of this equation. Hence we can define the reaction-rate densiry F(r,t) at any point in the system by
Furthermore the effort he expends in understanding the material in this chapter merely multiplying the neutron density N(r,t) by the interaction frequency uZ:
will provide him with a much deeper and more thorough understanding of-the
approximate methods we will develop in later chapters. expected rate at which
With this strong note of encouragement, let us now add a qualification. We have F (r, t)d3r IuZN (r, t) d3r= interactions are occurring (4-3)
attempted to develop these later approximations in a manner independent of this in d3r about r a t time t .
106 / THE ONESPEED DIFFUSION MODEL OF A NUCLEAR REACXOR NEUTRON TRANSPORT / 107
Notice that this "density" is defined with respect to both space and energy. One
can also generalize the concept of reaction rate density to include energy depen- [The term "angular" arises from the fact that n(r, E,fi,t) depends on* the velocity
dence as
+
spherical coordinate angles 0 and specifying the neutron direction a (see Figure
4-2).] This is the most general neutron density function we need to define since it
happens that one can derive an essentially exac: equation, the neutron transport
equation, for the angular neutron density n(r, E, a, t).
The product oN(r,t) arising in Eqs. (4-3) and (4-5) occurs very frequently in However before deriving this equation, it is useful to introduce several other
reactor theory, and therefore it is given a special name: definitions. We will first define the angular neutron flux in a manner similar to that
in which we earlier defined the neutron flux, simply by multiplying the angular
+(r, I ) = oN (r, I) r neutron flux [cm-'. sec- '1. (4-6) density by the neutron speed o:
Although it will certainly prove convenient to work with +(r,t) rather than N(r,t)
(since then one does not have to worry about including the neutron speed o in the
reaction rate densities), the tradition in nuclear engineering of referring to this A related concept is the angular current densiw, defined by
quantity as the neutron "flux" is very misleading. For +(r,t) is not at all like the
fluxes encountered in electromagnetic theory or heat conduction, since these latter
fluxes are wctor quantities, whereas +(r,t) is a scalar quantity. Actually the
"neutron current" J(r,t), which we shall introduce momentarily, corresponds more Notice that since d is a unit vector, the angular flux is actually nothing more than
closely to the conventional interpretation of a "flux." To avoid unnecessary the magnitude of the angular current density
confusion over this unfortunate convention, the student would probably do best at
this point to think of the neutron flux as simply a convenient mathematical variable
108 / THE ONESPEED DIFFUSION MODEL OF A NUCLEAR REACTOR
NEUTRON TRANSPORT / 109
Sometimes quantities such as N(r,t) and g(r,t) which do not depend on d are
referre! to as scalar *or total densities and fluxes, to distinguish them from
n(r, E, Sb, t ) and cp(r, E, St, 1). We find this nomenclature cumbersome and will avoid
it in our development.
Notice that if the angular density is independent of fi (i.e., it is isotropic) then we
find that Eq. (4-14) demands the presence of a 4 r normalization factor in the
angular density
The angular current density has a useful physical interpretation. Consider a small
area d A at a point r. [Here we will use the convention that dA=&,dA where 6, is and
the unit vector normal to the surface.] Then we can identify
/
FIGURE 4-3. Neutmpr incident w a dillerential
element of area dA . FIGURE 4 4 . Anisotmpies In the angular denslty n(r. E . G . ~ near
) a boundary or
a neutron source.
110 / THE ONESPEED DIFFUSION MODEL OF A NUCLEAR REACTOR NEUTRON TRANSPORT / 111
Notice that J(r,t) is actually what would be referred to as the "flux" in other fields 11. THE NEUTRON TRANSPORT EQUATION
of physics, since if we have a small area dA at a position r, then We will now derive an exact equation for the angular neutron density in a
system by simply balancing the various mechanisms by which neutrons can be
net rate at which neutrons pass
J(r, 1)-dA= (4-2 1) gained or lost from an arbitrary volume V within the system. That is, we will
through a surface area dA. consider mechanisms that will change the number of neutrons in this volume th5t
are characterized by a specific energy E and are traveling in a specific direction Q.
The units of both J(r.1) and +(r,t) are identical [~m-~.sec-'].However J is a It is convenient to use a bit of vector calculus here, but hopefully this will not
vector quantity that characterizes the net rate at which neutrons pass through a obscure the simple physics behind this. equation (which is just the mathematical
+
surface oriented in a given direction, whereas simply characterizes the total rate expression of a "count-the-neutrons" game).
at which neutrons pass through a unit area, regardless of orientation. Such an To this end, consider any old arbitrary volume-V. The number of neutrons in V
interpretation would suggest that J is a more convenient quantity for describing with energy E in dE and traveling in a direction S2 in dQ within this volume is just
neutron leakage or flow (e.g., through the surface of the reactor core), while + is
more suitable for characterizing neutron reaction rates in which the total number of
neutron interactions in a sample (e.g., a small foil) is of interest. Although the
angular flux and current density are very simply related, we will find that there is
no simple analogous relationship between J and +. These concepts may appear a (S!nce n(r, ~ , f i , i )is, a "density" in E and fi space, we must multiply it by dE and
bit confusing at first, but they will become more familiar after we have illustrated d P in order to get a number.) The time rate of change of this number, then, is
their application in both our further theoretical development and the problems at given by a balance relation
the end of the chapter.
A closely related concept is that of the partial current densities, J,(r, t ) which $ [J;(r,~,b,t)d'r 1d ~ d f i - g a i nin V-loss from V. (4-24)
correspond to the total rates at which neutrons flow through a unit area fromAleft
to right (J,) or right to left (J-). If we recall our earlier definition of j(r, E,Q,t), If we assume that the arbitrary volume V is chosen not to depend on time, we can
then it becomes apparent that bring the time differentiation inside the spatial integration
where 2n' is merely a convenient notation to indicate that the angular integration
is performed only over directions with components along the surface normal (2n+) We will now classify the various ways that neutrons can appear or disappear from
or in the opposite directi?n (2n-). For example, if we choose to define the polar V, and then we will try to write mathematical !expressions for each of these
coordinates that specify 51 along the normal to the surface, then in the integration mechanisms in terms of the angular density n(r, E, Q, 1).
for J + , + would range from 0 to 2n, while 0 would range only from 0 to n/2.
Gain mechanisms:
It is evident from this definition that
@Any neutron sources in V (e.g., fissions).
Neutrons streaming into VIthrough the surface S.
@Neu;rons of differtnt E', Q' suffering a scattering collision in V that changes
Hence J is sometimes referred to as the net current density, since it can be E', 51' into the E , n of interest.
constructed as the difference of the partial current densities.
then obviously
1
= [ ~ ~ 3 r v d ~ ~ n ( r , ~ . b . t ) d(4-35)
~db.
[This Lerm was really easy-we only needed to define a source density, Here we have noted that
s(r, E, a,0.1 v.vb=v6.v
@Loss due to collisions in V: The rate at which neutrons suffer collisions at a
point r is since d does not depend on r.
If we now combine all of these terms such that
Hence integrating this collision rate over the volume V , we find rate of change of number = 0 +Q +Q - - Q, (4-37)
of neutrons in V
then we find
@Gain due to neutrons scatiering into dE about E, dd about 6 from other
energies E' and directi%nsQ': lfi we recall from Chapter 2 that the probability
of scattering from E',Q' to E,Q is given in terms of the double-differeptial
scaitering cross section, then the rate at which neutrons scatter from E', n'to
E, Q is
However we now apply the fact that the volume V was quite arbitrarily chosen.
However we must consider contributions from-any E1,b'. Hence Hence the only way for the integral to vanish for any V is for its integrand to be
identically zero-that is,
LnY
" d 3 r j ( r ) = 0 * j(r)=O. (4-39)
This is known as the neutron transport equation. Several general features of the
equation should p e noted: First, it is a linear equation in the unknown iependent
variable n(r, E, Q, t ) with seven independent variables (r = x , y , z ; E; Q = O,+; t).
Since it contains both derivatives in space and time as well as integrals over angle
and energy, it is known as an "integrodifferential" equation.
However the presence of the derivatives suggest that we must also specify
appropriate initial and boundary conditions for the angular density. Since only a
single time derivative appears in the equation, we can simply choose the initial
condition to be the specification of the initial value of the angular density for all
positions, energies, and directions:
Initial condition: n (r,E, b,0 ) = n,(r, E, b),all r, E, 6 .
In this sense then, the system of algebraic equations we will arrive at for the
unknown components of f can be written as a marrix equation.
We must next replace the various operations in the original equation by their
discretized counterparts. For example, we would represent derivatives by finite
difrence formulas such as Actual function
where the wi are known as the quadrature weights. A thorough description of such
procedures can be found in any elementary textbook on numerical analysi~,'.~
Discrete representation
although frequently it is more userul to derive such numerical approximations
directly for the specific equation under investigation (as we will have occasion to
do in Chapter 5).
Such procedures lead eventually to a set of coupled algebraic equatiohs for the
components f , that can be solved on a digital computer. Frequently these discrete
values of the unknown f(x) provide an adequate representation. However occa-
sionally one wishes to reconstruct the original unknown f(x) for all values of x
from the discrete values inf. Then one must interpolate between the point valuesf,
at xi, for example by using polynomials. (See Figure 4-10.)
An alternative way to arrive at a discrete representation of an equation is to
write the unknown function as an expansion in a finite number of known functions
(frequently polynomials). If we call these expansion functions p,(x), then we would
write
Xo XI X, X3 ...
FIGURE 4-10. Discrete adinate rrpaeahIloa of a funclioo.
Hence once again we find that the function f (x) is represented by a vector
although in this case, the components of the vector are just the unknown expansion
coefficientsf,. Notice that if we can determine these expansion coefficients, then we
can easily reconstruct the unknown function j(x) by merely using Eq. (4-56).
Interpolation is not required as it is with the discrete ordinates approach. (Here we might recall Eq. (4-48) as an example in which such an expansion would
prove suitable.)
EXAMPLE: When the dependent variable ranges between - I and + 1, a very
convenient choice of expansion functions are the Legendre polynomials: There are a variety of techniques one can now use to obtain a set of algebraic
equations for the expansion coefficients from the original equation for f ( x ) . For
example, it is frequently possible to substitute the expansion Eq. (4-56) into the
original equation, multiply by each of the expansion functions p,(x), integrate over
This choice of expansion functions is frequently used to represent the angular the independent variable x, and then use various properties of thep,(x) (such as the
dependence of the neutron flux in one-dimensional problems in which p=cosO is property of orthogonality, which we will discuss later) to arrive at a set of algebraic
120 / THE ONE-SPEED DIFFUSZON MODEL OF A NUCLEAR REACTOR NEUTRON TRANSPORT / 121
equations for the f,. One can also use more elaborate schemes such as the calculus angular variables, one can use orthogonality to obtain a coupled set of equations
of oariarions or so-called weighted residual methods to arrive at the set of algebraic . for the expansion coefficients cp,(r,E,t). Since this set of equations is rather
equations. Since we will make only very limited use of such function expansions in
our elementary development of numerical analysis methods in this text, we will complicated when written out for general geometries, we will refer the interested
refer the interested reader to several detailed descriptions9~10
of these techniques for reader to other sources for the general form of the equation^.'^^
more information. In one-dimension, an expansion in spherical harmonics corresponds to an
expansion in Legendre polynomials, PI( p), where p = cos 0:
A. Discretization of the Angular Dependence
Let us first c~nsiderhow these techniques can be used to discretize the
direction variable Sl in the neutron transport equation. In the discrete ordinate
approach," we w o ~ l dfirst represent the independent variable by a discrete set of
:.,
directions or rays Q,,n = I,. N. In this case the general form of the equations for the expansion coefficients is
somewhat simpler and can be written as
We then represent functions of 5 by only their values at each of these mesh
directions:
1acp1 ( / + I )
--+--+-- acp1+1
I I + E m (x, E, r )
at (21+1) ax (zr+i) ax
where the w,, are appropriately chosen quadrature weights for the particular
numerical integration scheme used to handle the angular integrals. In this scheme,
the transport equation reduces to a coupled set of N equations of the form: This set of equations is known, naturally enough, as the P, equarions.3s4
In the particular case in which the expansion in spherical harmonics is truncated
after two terms, that is N = 1, the expansion for the angular flux takes the form:
the actual functional dependence of the neutron flux on the ind$pendent variable. Such sets of equations are known as the multigroup equations (in this case, the
The dependence of the angular flux on the neutron direction C4 is usually rather multigroup SN equations) and play a very important role in nuclear reactor analysis.
weak, hence a set of general functions such as the spherical harmonics will provide We will return in Chapter 7 to discuss their derivation in the case in which the
an adequate description. angular flux is treated within the diffusion approximation.
However, the neutron energy E spans an enormous range from lo-' eV up to
10' eV. The dependence of the neutron distribution on energy is determined by
quite different processes in different regions of energy. For example, at high
energies the neutron energy dependence is dominated by the fission spectrum. At C. Treatment of Space and Time
intermediate energies, neutron slowing down and resonance absorption are the
dominant processes, while at low energies, neutron thermalization is important. To
The final step is to discretize the space and time variables. Although this can
expand the neutron flux in a set of functions that adequately describe all of these
be done using function expansions (more specifically, expansions in the spatial or
processes is clearly hopeless. Indeed such function expansions are capable of
temporal modes of the system), such modal expansions are rarely used in general
describing neutron energy behavior only for a restricted range of neutron energies
studies (although they are occasionally used in certain types of specific calculations
or a specific reactor type. (An example of such an expansion known as energy or
in which the modes can be easily calculated). Instead one utilizes a direct discrete
spectrum synthesis is given in Chapter 13.)
ordinate treatment in most cases. First the spatial variables r = ( x , y , z ) are decom-
One must be careful even when applying the discrete ordinate approach to the posed into an appropriate spatial mesh. The various derivative terms are then
energy variable. The difficulties involved become quite apparent when the very replaced by finite difference equations defined on this mesh.
detailed dependence of the neutron cross sections on energy is recalled. It clearly
would be unthinkable simply to consider these cross sections tabulated at several
,,
Finally the time variable is broken into discrete time steps, say to,t t,, ... and the
corresponding time derivatives are replaced by suitable difference formulas. The
discrete points as an adequate representation of this detailed structure. detailed mesh structure and difference formulas one utilizes depends on the type of
Instead one first breaks up the neutron energy range into intervals or so-called problem one is investigating and will not be discussed here, since we will consider
energy groups: such topics in much greater detail in later chapters.
in order to calculate nuclear reaction rates and hence study the chain reaction (e.g.,
by calculating the multiplication factor k). So far everything is straightforward; but unfortunately the last term
evaluated in terms of Nr, E,t). In fact we find that a
a
cannot be
must be evaluated in terms
Surely we can formulate an equation for +(r, E , I ) by simply integrating the
transport equation over angle. Let's try and see what happens. That is, we w!l of the neutron current J by using Eq. (4-19):
integrate each term of the transport equation (4-43) over the direction variable a:
This is known as the neutron continuity equation, since it is just the mathematical
statement of neutron balance.
It is important to note that this equation contains two unknowns, Hr,E,t) and
We can simplify each of these terms somewhat by the straightforward manipula- J(r, E,?), unlike the neutron tr?nsport equation, which only contained one un-
tions indicated below: known, the angular flux p(r, E, P, I). Hence by removing the angular dependence
we have in the process introduced another unknown, J(r, E,t), and hence we now
have an insoluble problem (i.e., one equation in two unknowns). The moral of this
story is that you don't get something for nothing-that is, merely integrating out
the angular dependence doesn't remove the complexities of the angular variation. It
just shifts them by demanding that one obtain yet another equation relating
+(r, E, t) and J(r, E, 1).
It is impossible to express J(r, E,t) in terms of +(r, E,t) in a general and exact
manner. This is more apparent if we recall the definitions Eqs. (4-17) and (4-19):
0= i & ( r , ~,h,t=
) S(r, E, t). (4-75)
Here we have used our earlier expression Eq. (4-16) for the neutron flux 9 in terms
of the angular flux cp and also simply defined a sourceAter~ S(r, E, t). To evaluate
the inscattering term @ we first re:all ,!hat X,(E1-+E, W-tP) usually depends only
on the scattering angle cosine h=W.Q. This implies that
It is obvious that these two quantities are entirely different functions, although they
can boJh be expressed in terms of an angular integral of the angular flux
cp(r,E,P,t). Hence there is no reason why one would expect these functions to be
simply related.
where X,(E1-+E) is just the "single" differential scattering cross sectio? defin9 in Undaunted by our failure to find a simple equation for +(r, E, t), suppose we shift
Eq. (2-46). Hence we can interchange the order of integrations over D and 51' to our attention instead to developing an equation for the current density J(r, E, t). By
comparing the definitions in Eqs. (4-17) and (4-19) above, we are tempted to try
116 / THE ONESPEED DIFFUSION MODEL OF A NUCLEAR REACTOR NEUTRON TRANSPORT / 127
multiplying the transport equation by b a?d then integrating once again over so that @ finally becomes
angle. Actually since the direction variable 51 is a vector,
b=~,sin~cos~+~ysin~sin~+~zcos~,
A quick glance at the integral term confirms our fears; once again the streaming
term has kicked out yet another new unknown. To see this more clearly, we can
combine these results along with similar results for Qy and Q, to write Eq. (4-81) as
Each of these terms can be simplified in a manner similar to that used in deriving
an equation for the current density J:
the neutron continuity equation [Eq. (4-79))
However just as -with the neutron continuity equation [Eq. (4-79)], we find that
integrating over yields one equation but two unknowns, J(r, E,t) and
Now to handle the inscattering term @, write
[Here we t a r e taken the luxury of using a symbolic notation of writing two vectors
together, 5151. If this bothers you, just interpret this as a convenient notation for
taking each of the various combinations of components Q,Q,,Q,Q,, ...,a,&
vex! we do something a bit sneaky. Since b is a unit vector, we can write
separately to construct a quantity with nine components, II,,, n,, ...,n,,. Such
51'.CI1= 1. We will insert this into Eq. (4-85) so that we can rewrite it as:
quantities are referred to as tensors (or, in this case, dyadics), but we won't need to
get so formal here.] It should be evident that we can get a new equation for
s ( r , E,t) by multiplying the transport equation by fib and integrating, but this new
equation will contain yet another unknown,
Now we recall again- ht! X,(Ef-+E,h'+fi) depends only on the cosine of the
scattering angle h- a'. a. Thus we can write i p f i 6 b b o ( r , E, b. t).
4a I
. f i ) E1,bl.t) . (4-87)
@ = t i w d ~ * J _ d f i * [ j d~ ~ 5 1 ; ~ ~ ( E ~ + E , b 'Q;q(r, Hence all we are doing by multiplying by [fir
and integrating is generating an
infinite set of coupled equations (which we can't solve). [Incidently, jt should be
We will define apparent that the cylprit is the "streaming" or "leakage" term Q.Vq which
contains a factor of Q and hence generates the new unknowns in each equation.]
6. 11
d d i l x n ; ~ s ( ~ ~ - + ~ , b ~ . b )d=d b . r i q ( ~ * - + ~ , & . r i )
4n
The only way to cut off this chain of equations is to introduce an approximation.
We shall do this by assuming that the angular flux is only weakly dependent on
angle, and in so doing, we will generate the neutron diffusion equation. However it
is useful to first discuss several simplifications that can frequently be introduced
into the neutron transport equation.
128 / THE ONE-SPEED DIFFUSION MODEL OF A NUCLEAR REACTOR NEUTRON TRANSPORT / 129
B. Common Simplifications to the Neutron Transport Equation The assumption of isotropic neutron sources is usually not too restrictive since
most sources such as fission are indeed essentially isotropic. Unfortunately
1. THE ONE-SPEED APPROXIMATION
although neutron scattering is usually isotropic in the CM system, it is far from
It is frequently convenient to suppress the neutron energy dependence by isotropic in the LAB system, particularly for low mass number scatterers such as
assuming that one can characterize the neutrons by a single energy or speed. We hydrogen. Undeterred by such physical considerations, we will assume for the
will find in Chapter 7 that if one chooses the appropriate effective cross sections, moment that isotropic scattering is present. Then the one-speed transport equation
such a representation will in fact frequently yield a reasonable description of the simplifies still further to
reactor. However for now we will introduce the one-speed approximation in a rather
artificial manner by simply assuming that the neutron energy does not change in a
scattering collision. This can be inserted into the transport equation [Eq. (4-43)) in
a rather convenient manner by simply assuming a differential scattering cross
section of the form
However even this equation is extremely difficult to solve in general.
2,(E1-+E,Si'+fi) = X,(E,fi'+fi) 6(E1- E ) , (693) 3. OTHER SIMPLIFICATIONS
where 6(E1- E ) is the Dirac &-function defined by the property Thus far we have mutilated the energy and angular dependence of the
transport equation in the interest of mathematical expediency-and still have not
/ dx'f (xl)rS ( x - X I )-f ( x ) (4-94) arrived at anything we can hope to solve (at least analytically). So in frustration we
now turn our attention to the remaining time and spatial variables. First we will
for any sufficiently well-behaved function f(x). [See Appendix C for a more completely eliminate the time variable by agreeing to consider only steady-state
detailed discussion of animals such as the 8-function.] Using this definition, the transport problems. Then Eq. (4-99) simplifies to
inscattering term in Eq. (4-43) becomes
Next, we will assume that the system under study has uniform composition such
that the cross sections do not depend on position. Finally we will simplify the
Since all of the terms in the transport equation are now evaluated at the same system geometry, for example, by considering only planar or spherical symmetry.
energy, we may as well eliminate the explicit dependence on energy to write the In the case of planar symmetry we arrive at a rather simple-looking equation
one-speed neutron transport equation as
This equation is still far too complicated to solve (even using brute force numerical This equation can actually be solvedanalytically6-but only with rather sophisti-
techniques) in realistic geometries. So we'll introduce yet another simplification. cated mathematical techniques beyond the scope of this text. So even after a
number of rather questionable approximations, one arrives at an equation that can
2. ISOTROPIC SOURCES AND SCATTERING still only be solved with great difficulty.
After this rather pessimistic glance at the difficulties involved in solving the
One major simplification that can be introduced into the transport equation
arises when one assumes both isotropic neutron sources transport equation, let us remark that there is one very important class of transport
problems that can be solved exactly with only a minimal expenditure of effort-
those involving neutron transport in a purely absorbing medium.
4. NEUTRON TRANSPORT IN A PURELY ABSORBING MEDIUM
and isotropic scattering (in the LAB system) Frequently we are interested in neutron transport in a medium in which
scattering can be ignored. This might occur in a vacuum, for example (or more
realistically, a gas-filled region of a reactor). Or it might apply in a very highly
absorbing medium such as a fuel element or a control rod. In these cases, the
1.30 / THE ONE-SPEED DIFFUSION MODEL OF A NUCLEAR REACTOR NEUTRON TRANSPORT / I31
We have omitted the time dependence here since it is rarely of relevance when such
transport problems are of interest.
This equation can be solved exactly for any source distribution since it can be where 6(r-r') is just the three-dimensional version of the Dirac 8-function (see
converted into a simple first-order differential equation. Consider first the case of Appendix C) defined by
neutron transport in a vacuum in which X,=0:
Hence we find
However 6.0 is just the directional derivative in the direction 8. If we define a
variable R that measures distance along this direction (see Figure 4-1 1) then
This still looks a bit strange. [Actually it can be shown that p(r,6) vanishes unless
where R is measured in the - fi direction, and we find one is looking along the direction r, as one would expect from Figure 4-12.]
Suppose we compute instead the neutron flux itself:
Notice that this expression simply equates the neutron angular flux $1 position r in
direction d to the total nymber of source neutrons emitted in this $2 (obtained by
integrating back along - a).
Yet using the definition of 8(r-r') given by Eq. (4-108), we find that the flux
EXAMPLE: Consider an isotropic point source located at the origin (for con- resulting from an isotropic point source at the origin is just
venience, we will suppress the energy dependence for this example). The math-
that is, + falls off with distance as 1/rZ because of the ever increasing surface area
(4ar2) over which the So source neutrons/sec must be isotropically distributed.
,
FIGURE 4-11. Neutrw transport In a vacuum, I lsotmplenlly at the origin of an Infinhe medium.
132 / THE ONESPEED DIFFUSION MODEL OF A NUCLEAR REACTOR NEUTRON TRANSPORT / 133
However we can integrate this equation just as we did Eq. (4-103) to find
or canceling out the integrating factor from both sides FIGURE 4-13. CaMdiartes characterldngIdktwbed source s(r, ~,i).
We can finally use this result to synthesize the neutron flux resulting from an
arbitrary distribution of isotopic sources, S(r), in an infinite absorbing medium as
This soIution again has a v e j plausible interpretation when it is recognized that exp(- X,lr - r'l)
exp(- X,R) is just the attenuation that would occur between the source point and
the observation point r. Notice that this immediately reduces to our vacuum result
* 3 4nlr - f l 2 s (rl).
for the case in which Z,=O.
We can also obtain an exact solution for the situation in which 21, depends on
position. We need only use a slightly more complicated attenuation factor C. The One-Speed Diffusion Equation
We now turn our attention toward the development of an approximate
description of neutron transport more amenable to calculation than the neutron
where a(r,rl) is known as the optical thickness or optical depth of the media and is transport equation itself. To make life simple, we will first work within the
defined by one-speed approximation represented by Eq. (4-96). Let us first note the explicit
forms taken by the neutron conservation equation and the corresponding equation
for the current density J in the one-speed case:
EXAMPLE: Consider once again our point source, only this time assume that it
is imbedded at the origin of an infinitely large medium characterized by a uniform Here we have noted explicitly the simplifications that occur in the inscattering term
absorption cross section X,. Then repeating our earlier analysis using Eq. (4-1 15) when the one-speed approximation is introduced. More specifically,
yields
which is similar to our vacuum result, with the exception of an additional attenua- and
tion factor exp(-2,r) due to the absorption.
This very important result can be easily generalized to the situation in which the
source is located at an arbitrary point r' by merely shifting the coordinate system
origin to find But
134 / THE ONESPEED DIFFUSION MODEL OF A NUCLEAR REACTOR NEUTRON TRANSPORT / 135
where we have defined the average scattering angle cosine ji,, as and
As an aside, it should be noted that one can easily calculate & , for the case of However one can easily demonstraie (see Problem 4-14) that the integral of the
elastic scattering from stationary nuclei when s-wave scattering is present. For th$n product of any two components of gives
we know $at in the CM system, ac,(6c)Eos/4n. Hence if we use a c M ( ~ c ) d ~ c
= aL(B3dQL, we find
Hence
But recall
Of course Eqs. (4-132) and (4-135) are identical to our original definitions of the
1 + A cosBC flux and current earlier in Eqs. (4-17) and (4-19).
cos 6, =
We will now use the approximate form of the angular flux in Eq. (4-131) to
evaluate the second term in Eq. (4-122):
If we substitute this into Eq. (4-127) and perform the integration, we find the very
simple result
ji,,=-. 2 (4- 129) Next note that the integral of the product of any odd number of components of d
3A vanishes by symmetry:
Now that we have justified the forms of Eqs. (4-121) and (4-122) let us consider
how we might eliminate the annoying appearance of the third unknown,
j d b Bficp(r,fi,t). We will accomplish this by assuming that the angular flux is only
4 dq a ~ q
n
d -0 if i,n,or n is odd. (4-137)
weakly dependent on angle. To be more specific, we will expand the angular flux in
angle as If we use both Eqs. (4-134) and Eq. (4-137) we can evaluate
and neglect all terms of higher than linear order in fi. Actually a slightly different
notation for the unknown functions qo, qlx, ply, qII is useful. Write Eq. (4-130) as Hence by assuming that the angular flux depends only weakly on angle-more
specifically, that the angular flux is only linearly anisotropic-we have managed to
express the third unknown appearing in Eq. (4-122) in terms of the neutron flux
+(r,r). We have now achieved our goal of obtaining a closed set of two equations
1 3 for two unknowns, +(r, r) and J(r, 1):
=- +(r, t) + -J(r, 1).d . (4-131)
457 4n --
"
v at
+ Ve J + Z,,(r)+(r, I) = S (r, I), (4- 139)
Notice that we have labeled the unknown expansion coefficients as the flux and
current. That this notation is perfectly consistent can be seen by noting from Eq.
(4-13 1) that
I
1. ddq(r,d,t)-+(r,r);j;; J
4n
d d + ~ J ( r , t )=
. ( r , I), (4- 132) We have noted here that
Za(r) = Zt(9 - Zs(r)9
136 / THE ONESPEED DlFlWSlON MODEL OF A NUCLEAR REACTOR NEUTRON TRANSPORT / 137
and defined the macroscopic transport cross section quite frequently in other areas of physics where it is known as Fick's law. It is also
occasionally referred to as the diffusion approximation.
Before we consider the physical implications of this relationship, let us use it to
simplify the P I equations. If we substitute this into Eq. (4-139) we find
(We will comment on this definition in a moment.) These two equations are known
in nuclear reactor analysis as the P e uations (in the one-speed approximation)
I .9
since the approximation of linearly an~sotropic
angular dependence in Eq. (4-13 1)
in one-dimensional plane geometry is equivalent to expanding the angular flux in
Legendre polynomials in p =cosB and retaining only the 1=0 and I = I terms This very important equation is known,as the one-speed neutron d i f f i o n equation,
and it will play an extremely significant role in our further studies of nuclear
reactors. We will discuss its solution in considerable detail in Chapter 5, and we
will use it as the basis of a very simple but very useful model of nuclear reactor
hence the name PI approximation. Notice that this could be easily generalized to behavior.
obtain the P, approximation. Let us now return to consider the diffusion approximation [Eq. (4-148)] in more
In principle we could now use the P, equations to describe the distribution of detail. Notice that it implies that a spatial variation in the neutron flux (or density)
neutrons in a nuclear reactor. However it is customary to introduce two more will give rise to a current of neutrons flowing from regions of high to low density.
approximations in order to simplify thesz equations even further. First we will Physically this is understandable since the collision rate in high neutron density
assume that the neutron source term s(r,Q,t) is isotropic. This implies, of course, regions will be higher with the corresponding tendency for neutrons to scatter more
that the source term S,(r,t) vanishes in the equation for the current density. As we frequently away toward lower densities. The rate at which such diffusion occurs
mentioned earlier, this approximation is usually of reasonable validity in nuclear
reactor studies.
As our second approximation, we will assume that we can neglect the time
derivative o-' a J / a t in comparison with the remaining terms in Eq. (4-140). This
would imply, for example, that
I aiJl
- -<oZ,,
IJI at
that is, that the rate of time variation of the current density is much slower than the
collision frequency OX,.Since OX, is typically of order Id sec-I or larger, only an
extremely rapid time variation of the current would invalidate this assumption. We
will later find that such rapid changes are very rarely encountered in reactor
dynamics. Hence we are justified in rewriting Eq. (4-140) as
We can solve Eq. (4-145) for the neutron current density in terms of the neutron
flux
Hence we have found that in certain situations the neutron current density is
proportional to the spatial gradient of the flux. This very important relation arises
138 / THE ONESPEED DIFFUSION MODEL OF A NUCLEAR REACFOR NEUTRON TRANSPORT / 139
depends on the diffusion coefficient Which, in turn, is inversely proportional to the the energy-dependent P, equations
transport cross section Z,,.
It is convenient to introduce the concept of a transport mean free path
+ + ( E ) ,I ) J ~ ~ E ( E + E ) + (I ~ ,) + ( r ) (4-15))
at o
The transport mfp can be regarded as a corrected mfp accounting for anisotropies
in the scattering collision process. Since jb is almost always positive-that is, biased
in the direction of forward scarrering-the transport mfp A, will always be
somewhat larger than the actual mfp, AZ(X,)-'. This essentially accounts for the Continuing our analogy with the derivation of the one-speed diffusion equation, we
fact that neutrons experiencing forward scattering tend to be transported somewhat will again assume: (a) isotropic source S, =O (b) JJI-'aJJJ/at<oX,(r, E ) so that Eq.
further in a sequence of collisions than those being isotropically (or backward) (4-154) can be rewritten as
scattered.
Hence since D=&/3, we see that diffusion is also enhanced in a material with
pronounced forward scattering (e.g.. hydrogen), although this of course also de-
pends on the magnitude of the macroscopic cross sections.
It is important to keep in mind that the diffusion approximation is actually a However we can now see a problem that appears to prevent a straightforward
consequence of four different approximations: (a) the angular flux can be generalization of Fick's law to include energy dependence. For we cannot bring
adequately represented by only a linearly anisotropic angular dependence [Eq. J(r, E1,t) out of the scattering integral since it depends on the integration variable
(4-131)], (b) the one-speed approximation, (c) isotropic sources, and (d) the neutron E'.
current density changes slowly on a time scale compared to the mean collision time Of course, if we were allowed to assume isotropic scattering in the LAB system,
[Eq. (4-144)]. Actually only the first of these approximations is really crucial. The then & ( E '+E) =0 and we could find
remaining approximations can be relaxed provided we are willing to work with the
P I equations rather than the neutron diffusion equations (as one frequently is).
It is natural to ask when the angular flux is sufficiently weakly dependent on
angle so that the diffusion approximation is valid. More detailed studies of the
transport equation itself indicate that the assumption of weak angular dependence But the assumption of isotropic scattering is far too gross for most reactor
is violated in the following cases: (a) near boundaries or where material properties calculations.
change dramatically from point to point over distances comparable to a mean free We could proceed formally by merely defining an energy-dependent diffusion
path, (b)near localized sources, and (c) in strongly absorbing media. In fact strong coefficient
angular dependence can be associated with neutron fluxes having a strong spatial
variation. Usually if one is over several mean free paths from any sources or
boundaries in a weakly absorbing medium, the flux is slowly varying in space, and
diffusion theory is valid.
so that
If we now use this in the pair of Eqs. (4-79) and (4-91), we arrive immediately at ~ m d ~ ~ ~ s , ( ~ ~ + ~ ) ~ i ( ~ , ~ ~ ,(4-t I 60)
) = ~ z s ( ~ p i ( ~
140 / THE ONESPEED DIFFUSION MODEL OF A NUCZEAR REACTOR NEUTRON TRANSPORT / 141
Then we find a natural generalization of the diffusion coefficient: I . GENERAL MATHEMATICAL CONDITIONS ON THE FLUX
Although strictly speaking they are not boundary conditions, we should first
mention those mathematical properties that the function +(r,t) must exhibit in
order to represent a physically realizable neutron flux. For example, +(r,t) must be
Actually none of these derivations of an energy-dependent diffusion equation are a real function. Furthermore since both the neutron speed v and density N cannot
particularly satisfying because we have ignored the fact that in neutron transport be negative, we must require that +(r,t) be greater than or equal to zero. In most
processes, spatial transport, directional changes, and energy changes are intimately cases we can also require that +(r, t) be bounded. However we should add here that
mixed. We can only provide a more satisfactory derivation of Eq. (4-158) after we one occasionally encounters pathological models of physical neutron sources that
have discussed comparable approximations characterizing neutron energy transfer cause +(r,t) to diverge. An example would be the familiar point source. Another
in scattering collisions. Such a derivation must await our discussion of neutron situation is the so-called Milne problem, in which one studies the behavior of the
slowing down in Chapter 8. flux near a vacuum boundary fed by a source of infinite magnitude located at
For the present, we will simply assume that the generalization of Fick's law to infinity. There will also be certain symmetry conditions that we can place on +(r,t)
include energydependence is given by Eq. (4-158) with an energy dependent resulting from geometrical considerations, such as plane or axial symmetry. These
diffusion coefficient as defined by Eq. (4-161). If we now substitute this into Eq. conditions will become more understandable as we consider specific examples later
(4-l53), we arrive at the energy-dependent diffusion equation in Chapter 5.
2. BOUNDARIES AT INTERFACES
1 2-V .D (r, E)V+ +~ , ( rE, )+(*E, t)
v at
Consider next an interface between two regions of differing cross sections.
Now clearly the correct transport boundary condition is that
E )qi(r, E ',t)
= I w d E 1z~(E'+ + S(r, E, t).
cp,(r,,b, t) = cp2(rs,8$1)for all 8 ,
This equation plays a very important role in nuclear reactor analysis since it is
frequently taken as the starting point for the derivation of the multlgroup diffusion where cp, is the angular flux in region 1, while cp, is the angular flux in region 2.
equations. These latter equations represent the fundamental tool used in modern This condition ensures conservation of neutrons across the boundary and can
nuclear reactor analysis. easily be derived directly from the transport equation.
Unfortunately we cannot satisfy this boundary condition exactly using diffusion
E. Diffusion Theory Boundary Conditions theory. At best, we can only ensure that angular moments of Eq. (4-166) are
satisfied. And since diffusion theory yields only the first two moments of the
Since the neutron diffusion equation has derivatives in both space and time, it angular flux, +(r,t) and J(r,t), at best we can demand
is apparent that one must assign suitable boundary and initial conditions to
complete the specification of any particular problem. Since the diffusion equation
itself is only an approximation to the more exact transport equation, we might
suspect that we can use the transport theory boundary conditions as a guide in our
I,.t p I r , , t = J b ( r , , ) + l ( s ) = s
4n
(4-167)
We can obtain the appropriate initial condition for the diffusion equation by
merely integrating the transport condition over angle to obtain:
and Once again diffusion theory will only be able to approximate this boundary
condition. Notice in particular that the boundary condition is given only over half
of the range of solid angle (corresponding to incoming neutrons). Hence suppose
we again seek to satisfy the transport boundary condition in an "integral" sense by
Hence the interface diffusion theory boundary conditions are simply those corre- demanding
sponding to continuity of flux and current density across the interface:
where 6, is the unit normal to the surface. Hence our diffusion theory approximation to the transport boundary condition
Eq. (4-171) is just
3. VACUUM BOUNDARIES
Recall that our transport theory boundary condition was merely a mathemati-
cal statement that there could be no incoming neutrons at a free or vacuum
boundary For convenience consider this boundary condition applied to a one-dimensional
geometry with the boundary at x = x,
cp(r,,b,t)==o for djl.ds<o,
all r, on S.
Notice that this relation implies that if we "extrapolated" the flux linearly beyond
the boundary, it would vanish at a point
For this reason, one frequently replaces the vacuum boundary condition
J-(x,)=O, (4- 177)
REFERENCES
4-1 We have defined the angular neutron density n ( r , ~ , b , t )in terms of the neutron
energy E and the direction of motion 6 , but one could as well define an angular
More complicated extrapolation length formulas can be derived for curved density that depends instead on the neutron velocity v, n(r,v,t). Calculate the
boundaries," although they are rarely necessary since Eq. (4-180) suffices unless relationship between these two dependent variables.
the radius of curvature of the boundary is comparable to a mean free path. 4-2 Two thermal neutron beams are injected from opposite directions into a thin sample
It should be remembered that the true flux does not vanish even outside the of 235U.At a given point in the sample, the beam intensities are 1012neutrons/cm2.
boundary. The diffusion theory flux is a poor representation of the true flux near sec from the left and 2 x 10" neutrons/cm2.sec from the right. Compute: (a) the
the boundary (as we saw earlier, diffusion theory is not valid near a boundary). neutron flux and current density at this point and (b) the fission reaction rate density
The boundary conditions we have derived are intended to yield the proper flux at this point.
only in the interior of the reactor, that is, several mean free paths away from the 4-3 Suppose that the angular neutron density is given by
reactor boundary.
These boundary conditions complete our'description of neutron transport within where 0 is the angle between d and the z-axis. If A is the area perpendicular to the
the diffusion approximation. The neutron diffusion equation will play a very z-axis, then what is the number of neutrons passing through the area A per second:
fundamental role in our development of nuclear reactor analysis methods. We will (a) per unit solid angle at an angle of 45" with the z-axis, (b) from the negative z to
begin our study of nuclear reactor behavior using one-speed diffusion theory, since the positive z direction, (c) net, and (d) total?
while this description has only limited quantitative validity, it does allow us to 4 4 In a spherical thermal reactor of radius R , it is found that the angular neutron flux
can be roughly described by
illustrate rather easily the principal concepts of nuclear reactor theory, as well as to
develop the mathematical techniques used in more sophisticated models. With this
background, we will then develop and apply the principal tool of the nuclear
reactor designer, multigroup diffusion theory. Compute the total number of neutrons in the reactor.
146 / THE ONE-SPEED DIFFUSION MODEL OF A NUCLEAR REACTOR NEUTRON TRANSPORT / 147
4-5 Demonstrate that in an isotropic flux, the partial current density in any direction is 4-12 Use Simpson's rule to write a numerical quadrature formula for the angular integral
given by J + = +/4. j ? ~dpv(x,p)for N equal mesh intervals.
4-6 One of the principal assumptions made in the derivation of the neutron transport
4-13 Develop the multigroup form of the transport equation as follows: First break the
equation is that there are no external forces acting on the neutrons, that is, that the energy range 0 < E < 10 MeV into G intervals or groups. Now integrate each of the
neutrons stream freely between collisions. Suppose we were to relax this assumption terms in the transport equation, Eq. (4-43), over the energies in a given group, say
and consider such a force acting on the neutrons (say, a constant gravitational force E, < E < E,- ,. (Remember that the group indexing runs backwards such that O = E,
in the z-direction). Then we might expect additional terms to appear that involve this < EG- I < - .. E8 < E8- I < .- . El < Eo= 10 MeV.) Now by defining the groupJuxes as
force. Derive this more general transport equation. the integral of the flux over each group, and the cross sections characterizing each
4-7 Explain briefly whether or not the transport equation as we have derived it adequately group as in Eq. (4-69), determine the set of G equations representing the transport
describes the spatial and angular distribution of neutrons: (a) in a flux of the order of equation.
100/cm'.sec, (b) scattering in a single crystal, (c) passing through a thin pure 4-14 By writing out the components of the direction unit vector db in polar coordinatu,
absorber, and (d) originating within a very intense nuclear explosion in which the demonstrate explicitly that
explosion debris are moving outward with very high velocities.
4-8 Develop the particular form of the transport equation in spherical and cylindrical
geqmetriu. To simplify this calculation, utilize the one-speed form of the transport
equation (4-100) in which time dependence has been ignored.
4-9 An interesting model of neutron transport involves transport in a one-dimensional
rod. That is, the neutrons can move only to the left [say, described by an angular
density n-(x, E,t)] or to the right [n+(x,E,r)]. One need only consider forward
i w d dQi -0 and 14r d d QIR, 3
I-
4n
3
0
i#j
scattering events described by Z,+(Ef-+E) or backward scattering events described by Verify Eq. (4-138) using the identity in Problem 4-15.
Z;(E'-+E). Perform the following: (a) derive the transport equations for n + and n-, Explicitly demonstrate by integration that p,,=2/3A for elastic scattering from
(b) make the one-speed approximation in this set of equations, that is, stationary nuclei when such scattering is assumed to be isotropic in the CM system.
Consider an isotropic point source emitting So monoenergetic neutrons per second in
z,(E)-~(E)+z:(E)+z;(E)=z,+I:+~;, an infinite medium. Assume that the medium is characterized by an absorption cross
section Z,, but only by negligible scattering. Determine the rate at which neutron
absorptions occur per unit volume at any point in the medium.
Compute the neutron flux resulting from a plane source emitting neutrons isotropi-
and (c) describe the boundary and initial conditions necessary to complete the cally at the origin of an infinite absorbing medium. Hint: Just represent the plane
specification of the problem if the rod is characterized by a length L. source as a superposition of point sources.
4-10 Consider the following differential equation: In a laser-induced thermonuclear fusion reaction, a tiny pellet is imploded to super
high densities such that it ignites in a thermonuclear bum. In such a reaction some
d2f 10" 14 MeV neutrons will be emitted essentially instantaneously (within 10-I' sec).
- +xY(x)-2x(4- x),
Compute the neutron flux at a distance of 1 m from the reaction as a function of time,
dx2
assuming that the chamber in which the reaction occurs is evacuated.
where j ( x ) is defined on the interval 0 < x <4. Discretize this equation by first Compute the thermal neutron diffusion coefficients characterizing water, graphite,
breaking up the independent variable range into four segments of equal length. Next and natural uranium. Then compute the extrapolation length zo characterizing these
use a finite-difference approximation to the d2f/dx1 term to rewrite the differential materials.
equation as a set of algebraic equations for the discretized unknown f(xl)=f,. For Assuming that the diffusion approximation [Eq. (4-131)] is valid, compute the partial
convenience, assume boundary conditions such that f(0) = 0- f (4). Solve this set of current densities in the z direction defined by Eq. (4-22). Use the one-speed
equations for the/, and then plot this solution against x using straight-line interpola- approximation.
tion. Consider the time-dependent one-speed P I equations assuming isotropic sources and
4-1 1 Consider the steady-state one-speed transport equation assuming isotropic scattering plane symmetry. Eliminate the current density J(x,t) to obtain one equation for the
and sources in a one-dimensional plane geometry neutron flux $~(x,t). Compare this equation with the one-speed neutron diffusion
equation and indicate what differences you might expect in solutions to the two
equations.
4-26 Try to construct solutions to the one-speed transport equation in an infinite sourceless
medium:
Expand the solution to this equation in the first two Legendre polynomials
Substitute this expansion into the transport equation, multiply by P,(p) and Pl(p)
respectively, and integrate over p to obtain a set of equations for the unknown by seeking solutions of the form p ( x , p ) = ~ ( p )exp(- x/v) where v and ~ ( pare
) to be
expansion coefficients cpdx) and cpl(x). (These are just the P I equations.) determined.
148 / THE ONE-SPEED DIFFUSION MODEL OF A NUCLEAR REACTOR
4-27 Consider the time-independent one-speed transport equation under the assumption of
isotropic sources and scattering
The One-Speed
Diffusion Theory Model
By regarding the right-hand side of this equation as an effective source, use the result
[Eq. (4-120)] to derive an integral quation for the neutron flux Hr).
A. Derivation of the Diffusion Equation Hence the time rate of change of the number of neutrons in V must be just
We will suppress the neutron energy dependence by assuming that all of the
neutrons can be characterized by a single kinetic energy. Such a one-speed (or
one-group) approximation greatly simplifies the mathematical study of nuclear
reactor behavior. Of course such an approximation is also highly suspect, particu-
= Production in V - absorption in V
larly in light of the fact that neutron energies typically encountered in a reactor
span a range from to 10' eV, and neutron cross sections depend sensitively on -net leakage from V . (5-2)
energy over most of this range. We will later be able to show, surprisingly enough,
that if one regards the one-speed approximation as an average description and We can easily write down mathematical expressions for the gain and loss terms. If
chooses the appropriate average cross sections, then in fact the one-speed model we define a neutron source density S(r,t), then
can actually be used to obtain a quantitative description of a nuclear reactor.
As yet a further justification for our exhaustive study of the one-speed approxi- Production in V - J:'r s (r,t).
mation, we would remark that most energy-dependent theories (e.g., multigroup
diffusion theory) are solved by performing a sequence of one-speed calculations for
each successive energy group. Hence the methods we develop for analyzing our Since the absorption rate density at any point in V is just X,(r) +(r,r), it is obvious
one-speed model will later be extended directly to more sophisticated descriptions. that the total rate of neutron loss due to absorption in V is just
We will characterize the neutron distribution in the reactor by the neutron
density N(r,r) which gives the number of neutrons per unit volume at a position r Absorption in V = p rXa (r)@(r,t).
at time t. Actually we will find it more convenient to work with the neutron flux,
+(r,t)==vN(r,r), since then we can compute the rate at which various types of The term describing neutron leakage out of or into V is a bit more difficult. If
neutron-nuclear reactions occur per unit volume by merely multiplying the flux by J(r,t) is the neutron current density, then the net rate at which neutrons pass out
the corresponding macroscopic cross section. For example, the rate at which fission through a small surface element d S at position r, is J(r,,r).dS. Hence the total net
reactions occur per cm3 at a point r would be given by Ll,(r)+(r,t). leakage through the surface of V is just
We will derive an equation for the neutron flux by merely writing down a
mathematical statement of the fact that the time rate of change of the number of Net leakage from V = I d s . J(r,r). (5-5)
S
neutrons in a given volume must be simply the difference between the rate at which
neutrons are produced in the volume and the rate at which they are lost from the
Now we could combine all of these terms back into Eq. (5-2) as they stand, but first
volume due to absorption or leakage.
it is convenient to convert Eq. (5-5) into a volume integral similar to the other
terms. The common way to convert such surface integrals into volume integrals is
to use Gauss's theorem to write
If we now substitute each of these mathematical expressions into Eq. (5-2), we find
But recall that we chose our volume .V to be arbitrary. That is, Eq. (5-7) must hold equation
for any volume V that we would care to examine. However the only way this can
occur is if the integrand itself were to vanish. Hence we find we must require
This equation will form the basis of much of our further development of nuclear
reactor theory.
Of course this equation is still quite exact (aside from the one-speed approxima-
tion), but it is also quite formal since it contains two unknowns, +(r,t) and J(r,t). B. Initial and Boundary Conditions
To proceed further, we need yet another relationship between these two variables.
Unfortunately there is no exact relationship between +(r,t) and J(r,t). One must We must augment this equation with suitable initial and boundary conditions.
resort to an approximate relation. Now it is well known from other fields of physics Although these conditions have been developed in a more rigorous fashion from
such as gaseous diffusion that the current density is approximately proportional to neutron transport theory in Chapter 4, we will remotivate them here using plausible
the negative spatial gradient of the density, or in our case, the flux. That is, physical arguments.
particles will tend to flow from regions of high to low density at a rate proportional The appropriate initial condition involves specifying the neutron flux +(r,t) for
to the negative density gradient. Stated mathematically, one finds that all positions r at the initial time, say t =O:
Initial condition: +(r,O) = +,(r), all r. (5- 12)
The boundary conditions are a bit more complicated and depend on the type of
where the' constant of proportionality D (r) is known as the d~$usion coefficient, physical system we are studying. The principal types of boundary conditions we
while Eq. (5-9) is referred to as Fick's Imu. Of course to postulate such a will utilize include the following:
relationship between current density and flux implies nothing about its range of
validity. Indeed we do not even know what the diffusion coefficient D is. This 1. VACUUM BOUNDARY
situation is very common in macroscopic descriptions of physics, and relationships
such as Eq. (5-9) are usually referred to as "transport laws" (not to be confused At the outside boundary of a reactor, one would like to construct a boundary
with the transport equation-we are talking the language of the physicist now) condition corresponding to the fact that no neutrons can enter the reactor through
while the proportionality coefficients are known as "transport coefficients." Ex- this surface from outside. Implicit in this fact is the assumption that the reactor is
amples are numerous and include Fourier's law of thermal conduction (thermal surrounded by an infinitely large vacuous region. Of course no reactor is
surrounded by a vacuum, but rather by air, concrete, and a host of other materials.
conductivity), Stokes' law of viscosity (shear viscosity), Ohm's law (electrical resis-
It is frequently convenient to assume that the reflection of neutrons back into the
tivity), to name only a few. In all cases, one is forced to go to a microscopic
description in order to evaluate the transport coefficient and examine the range of reactor from such materials is negligible so that nonreentrant boundary conditions
validity of the transport law. apply.
However this is of course exactly what we did by deriving the diffusion approxi- There is only one problem; diffusion theory is incapable of exactly representing
mation [Eq. (5-9)) from the neutron transport equation in Chapter 4. There we a nonreentrant boundary condition. The closest one can come would be to demand
found that that the inwardly directed partial current
where & is the average cosine of the scattering angle in a neutron scattering
collision. Furthermore, we found that Eq. (5-9) was valid provided it was used to vanish on the boundary (clearly an approximation, since this expression for J -
describe the neutron flux several mfp away from the boundaries or isolated already is approximate). Actually we really shouldn't worry much about a consis-
sources, the medium was only weakly absorbing, and provided the neutron current tent free surface boundary condition within the diffusion approximation, for we
was changing slowly on a time scale comparable to the mean time between have already indicated that diffusion theory is not valid near the boundary anyway.
neutron-nuclei collisions. It is important to keep these limitations in mind as we It can only be expected to hold several mfp inside the boundary.
apply this approximation in nuclear reactor analysis. Hence what we should really look for is a "fudged-up" boundary condition
Henceforth we will accept the diffusion approximation [Eq. (5-9)] as providing a which, although it may have little physical relevance at the boundary, does in fact
valid expression for the neutron current density in terms of the neutron flux. If we yield the correct neutron flux deep within the reactor where diffusion theory is
substitute Eq. (5-9) into Eq. (5-8). we arrive at the one-speed neutron diffurion valid. More detailed transport theory studies indicate that the proper boundary
I54 / THE ONE-SPEED DIFFUSION MODEL OF A NUCLEAR REACTOR THE ONESPEED DIFFUSION THEORY MODEL / 15s
condition to choose is one in which the diffusion theory flux g(r,t) vanishes at a occasionally encounter a situation in which the neutron flux becomes singular at a
distance localized source (e.g., a point source). However since such sources are mathemati-
cal idealizations, this singular behavior doesn't bother us, and in general we will
demand that the flux be finite away from such sources. This condition is particu-
outside of the actual boundary of the reactor. This extrapolafed boundary is usually larly useful in geometries in which certain dimensions are infinite.
denoted symbolically with a tilde. For instance, if the physical surface is at r,, then We will also occasionally be able to use symmetry properties to discard physi-
the flux will be taken to vanish on the extrapolated boundary qb(?,,r)=O. cally irrelevant solutions of the diffusion equation. For example, in one-
A side comment here is appropriate. For most reactor materials, A,, is quite dimensional slab geometries, we can choose the coordinate origin at the centerline
small, usually being of the order of several centimeters or less. When it is of the slab, and then use symmetry to eliminate solutions with odd parity [i.e.,
recognized that most reactor cores are quite large (several meters in diameter), then $4 - x) = - cP(x)l.
it is understandable that one frequently ignores the extrapolation length zo and Other types of boundary conditions are encountered in practice. A very common
simply assumes the flux vanishes on the true boundary. problem in reactor calculations involves the determination of the flux in a small
Furthermore few realistic reactor geometries are surrounded by a free surface. subregion of the reactor fuel lattice, a so-called unit cell repeated throughout the
lattice. For example, such a cell might contain a single fuel rod surrounded by
Rather they are surrounded by coolant flow channels or plenums, structural
coolant (a fuel cell) or several fuel assemblies along with a control element (a
materials, thermal shields or such. Hence while the concept of a free surface is
control cell). Since these unit cells are repeated in a regular fashion throughout the
useful pedagogically for painting a picture of an idealized reactor geometry
core lattice, one can argue that there should be no net transfer of neutrons between
surrounded by a vacuum of infinite extent, it is rarely employed in modern nuclear cells, that is, that the neutron current density J(r) vanish on the cell boundaries.
reactor analysis. This is an example of a boundary condition on the current. In such cell calcula-
2. INTERFACES (MATERIAL DISCONTINUITIES) tions it is also frequently necessary to obtain a diffusion theory solution in the
vicinity of a strong absorber (e.g., a fuel rod or control element). The appropriate
The structure of a nuclear reactor core is extremely complex, containing boundary condition at the interface between the diffusing medium and the ab-
regions of fuel, structural material, coolant, control elements, and so on. Hence sorber is handled much like that characterizing a free surface. That is, one uses a
while one rarely encounters situations in which the material cross sections Z(r) transport-corrected boundary condition on the flux or current to yield the proper
depend continuously on position, one frequently is faced with what might be diffusion theory solution at a distance of several mfp from the interface. We will
termed "sectionally uniform" cross sections that change abruptly across an inter-
face separating regions of differing material composition. Our usual procedure in
treating such discontinuities in material properties will be to solve the diffusion ,Equivalent fuel cell
equation in each region separately and then attempt to match these solutions at the
interface using appropriate boundary conditions. Moderator
/ Fuel slemnt /Control rod /Fuel anerblv
Once again the diffusion equation is not strictly valid within several mfp of the Fuel element
interface. However in this case one can argue that conservation of neutron
transport across the interface demands continuity of both the normal component of
the neutron current density J(r,r) and the neutron flux &(r,r).
This condition is occasionally modified to account fb; tl;e physically fictitious
but mathematically expedient convenience of including an infinitesimally thin
neutron absorber or source at the interface. Then while the neutron flux i s still
continuous across the interface, the normal component of the current experiences a
jump:
where $ is the interface surface normal, while S would represent a source term (if
positive) or an absorption term (if negative). (See Figure 4-15.)
give explicit examples of such boundary conditions when we consider cell calcula- We will frequently consider situations in which the flux is not a function of time.
tions in Chapter 10. Then Eq. (5-19) becomes
We should reemphasize that the application of neutron diffusion theory in
reactor analysis is successful in large part because the diffusion equation and its +
- D V2+(r) Z,+(r) = S (r).
boundary conditions have been modified using more exact transport theory correc-
tions. For example, D=[3(X,-&Z,)]-' contains a correction to account for This equation is known as the Helmholtz equation and is also a very familiar beast
anisotropic scattering. Furthermore the boundary conditions on the flux or current in mathematical physics. It is useful to divide by - D to rewrite Eq. (5-20) as
at a free surface or adjacent to a highly absorbing region contain transport
corrections to yield the proper neutron flux deeper within the diffusing region.
Such transport corrections frequently yield diffusion theory estimates that are far
more accurate than one would normally expect, especially when we recall the
rather strong approximations required to derive the neutron diffusion equation. where we have defined the neutron diffusion length L
slightly weird source. Note that if we restrict x#O, the source term disappears from
Eq. (5-23):
Our approach will be to solve this homogeneous equation for x # 0, and then use a Hence our solution is
boundary condition at x = O to "fix up" these solutions.
We could obtain this boundary condition directly by integrating Eq. (5-23) from
x = 0- c to x =O+ c across the source plane and then taking the limit as € 4 0 to
find and by symmetry we can infer
[See Problem 5-21. However we might also merely note that this is just a special
case of the more general interface boundary condition [Eq. (5-15)]. If we use the Hence the neutron flux falls off exponentially as one moves away from the source
symmetry of the geometry to assert that J,(O+)= -J,(O-)=J(O), then our plane with a characteristic decay length of L. As one might expect, the larger L
boundary condition at the source plane becomes just (i.e., the smaller Z,), the less the neutron flux is attenuated as we move into the
medium. Notice also that the magnitude of the flux is proportional to the source.
so
lim J ( x ) = - That is, doubling the source strength will double the neutron flux ' ( x ) at any
x+o+ 2 position x , but this should have been anticipated since the neutron diffusion
This source boundary condition makes sense physically, since it merely says that equation is linear and hence the principle of superposition holds.
the net neutron current at the origin on either side must be just half of the total The reader should be reminded that, while this solution may provide a reason-
source strength. able description of the neutron flux in a medium (provided it is not too highly
We are not through with boundary conditions yet. Since we have a second-order absorbing), it is certainly not valid within several mfp of the source plane itself.
derivative, d 2 / d x 2 , we need another boundary condition. We will use the boundary Indeed more accurate transport theory studies6*'indicate that the neutron flux does
condition of finite flux as x+oo. not look at all exponential near the source. In fact there are additional components
Hence the mathematical problem to be solved is to the solution. Fortunately if absorption is not too strong, these transport theory
"transients" rapidly diminish as one moves several mfp away from the source
plane, and the simple exponential behavior predicted by diffusion theory is found,
with one mild modification. Transport theory predicts that the "relaxation length"
with boundary conditions: characterizing the exponential decay is not L = ( D / z ~ ' / ~but
, rather is given by
the root of the transcendental expression
lim - D - d'= - so
(a) X+O+ dx 2
(b) J& '(x) <
Fortunately if Z,<Z, (as it must be for diffusion theory to be valid), one can
160 / THE ONE-SPEED DIFFUSION MODEL OF A NUCLEAR REACTOR THE ONESPEED DlFrmSlON THEORY MODEL / 161
For example, in graphite X,= .385 cm-I while 2,s .00032. Hence X,/X,= 8.3 x Hence the neutron-diffusion length L has the interesting physical interpretation as
lo-', which implies that the correction to L given by higher terms in this expansion being I / Gof the root mean square (rms) distance to absorption
(that is, by transport theory) is only about 0.03% in this material.
Applying the boundary conditions, we find that (b) implies that we choose B = 0 ,
162 / THE ONESPEED DIFFUSION MODEL OF A NUCLEAR REACTOR THE ONESPEED DIFFUSION THEORY MODEL / 163
where the integral is taken over all space. Then boundary condition (a) implies
The cylindrical geometry problem of a line source at the origin of an infinite
medium can be worked out in a very similar way. However to do so here would
"deprive" the reader of an opportunity to try his own hand at such diffusion theory
problems. Hence we have left the line source as an exercise in the problem set at
the end of the chapter. We will instead turn our attention to problems in finite Our final solution is therefore
geometries.
Here we have replaced the boundary condition at infinity by the vacuum boundary
condition-in this case, using an extrapolated boundary Z / 2 = a / 2 + 2,. If we again where L, ={% is the diffusion length characterizing region @. Similarly in
seek a general solution of the form of Eq. (5-28), then applying the boundary region @
condition (b) implies
Now we will need several boundary conditions. We can use our earlier conditions been lost to the vacuum. Such materials used to reduce neutron leakage are known
IEq. (5-2711 as reflectors. Any material with a large scattering cross section and low absorption
cross section would make a suitable neutron reflector. For example, the water
lim J,(x) = 7
so channels surrounding LWR cores act as reflectors. In the HTGR, graphite blocks
(a) x+o+ are added to the top and bottom of the reactor core to serve as neutron reflectors.
(b) +2(~)<w as x-oo. A concept very closely related to that of neutron reflectors is the rejection
coefficient or albedo, defined as the ratio between the current out of the reflecting
In addition we will use the interface conditions region to the current into the reflecting region:
JO",
a=-.
Jin
To make this concept more precise, suppose we want to attach a reflecting slab of
thickness a to a reactor core (or perhaps a medium with a neutron source in it such
Using our earlier work as a guide, we can seek general solutions as the slab geometry we have just considered). If we are given a current density
Jim= J + entering the reflecting slab surface, which we locate at x=O (see Figure
5-7) for convenience, then we can solve the diffusion equation characterizing the
+,(x) = A , cosh -
X
+
B, sinh 2 in region @ reflecting region:
- Ll
+2(x) ~ , e x p ( -
L1
%) + B 2 e x p ( ~ ) in region @
and apply the boundary conditions to find (after a bit of algebra)
subject to boundary conditions J+(O)= Jin,+(Z)=O. We can then solve for the flux
+(x) in the reflecting slab and use this solution to compute the albedo a as (see
Problem 5-1 I):
It is of interest to plot the albedo for the slab reflector versus slab thickness as
shown in Figure 5-8. For thin reflectors, very few of the neutrons are reflected and
hence the albedo is small. As the reflector becomes very thick, the albedo
We have sketched the form of this solution in Figure 5-6. Several features of this
solution are of some interest. Note that while the neutron flux is continuous across
the interface, the derivative of the flux is not. This later discontinuity is, of course,
a consequence of the fact that the diffusion coefficients in the two regions differ;
hence to obtain continuity of current J , we must allow a jump in d+/dx across the
interface.
We have compared the solution for this problem with that obtained earlier for
the slab surrounded by a vacuum. It should be noted that the flux in the central
region falls off somewhat more slowly when the vacuum is replaced by a diffusing
material. This can be readily understood by noting that the material surrounding
the slab will tend to scatter neutrons back into the slab that would have otherwise
166 / THE ONESPEED DIFFUSION MODEL OF A NUCLEAR REACTOR THE ONESPEED DlFFUSlON THEORY MODEL / 167
Suppose this source was located at the point r' instead. Then the flux could be
found by a simple coordinate translation as
Y a
FIGURE 5-8. I%e alkdo f a r llalte shb plotted
versos s*b tblckae~.
approaches an asymptotic limit dependent only upon the material properties D and Next suppose we have several point sources at positions r;, each of strength Si.
L: Then we can use the fact that the diffusion equation is linear to invoke the
principle of superposition and write
Ir - rjl
sexd- -\
For example, the albedos characterizing infinitely thick reflectors of graphite, H,O,
and D,O are 0.93. 0.82, and 0.97 respectively.
The albedo can be used to replace the detailed solution in the reflecting region
by an equivalent boundary condition at the edge of the core using Eq. (5-45). That
is, if we use our earlier definitions of the partial current densities J , , then the
effective boundary condition on the flux in the region x <O is just
A n g l e point source
Used in this manner, the albedo becomes a very useful device for obtaining
boundary conditions for reactor core calculations.
One can continue this game of solving the diffusion equation in various one-
dimensional geometries indefinitely. As we mentioned earlier, it is simply an
exercise in ordinary differential equations. A variety of two- and three-dimensional
problems can also be studied. However since these latter problems involve partial
/ & e r a ~ point sources
differential equations, we will defer their treatment until after we have introduced
the separation of variables approach for solving such equations in Section 5-111-C.
If we really want to get masochistic, we can remove some of the symmetry in our
earlier problems so that the original partial differential equation (5-20) would have
to be solved directly-for example, a point source set off-center in a sphere.
However there is very little in the way of new physics to be learned from such
exercises. Hence we will bypass further examples in favor of moving directly to
more general problems. In particular, we will study how our previous solutions for
plane and point sources can be used to determine the neutron flux resulting from
an arbitrary distribution of neutron sources. /istributed source FIGURE 5-9. Supuposltloaof several point sounos.
168 / THE ONESPEED DIFFUSION MODEL OF A NUCLEAR REACTOR
THE ONESPEED DIFFUSION THEORY MODEL / 169
Once again we will allow the reader the privilege of deriving the line source
diffusion kernel. Actually both the plane source and line source kernels could have
This is frequently rewritten as been superimposed from the point source diffusion kernel.
Here K(x,x') would be the known kernel, while f(x) might either be known [as the
source density S(r)] or unknown as the flux cp in the inscattering term of the
transport equation
(We will leave it as understood that one must also apply suitable boundary
conditions.) There are a variety of techniques available to solve such inhomo-
geneous problems. The approach we have been using thus far is known as the
Green's function m e t h ~ d : ~
where x ~ ( E ' +E, d1+6)is known as the scattering kernel.] (a) GREEN'S FUNCTION METHODS
As a second example of such kernels, consider the flux resulting from a plane In this technique we first construct the solution to
source at the origin
M+8(x) = 6 (x - x'). (5-64)
Then if we call +g(x)= G(x,xl) the "Green's function" for the operator M, we find
that the general solution to Eq. (5-64) is just
If this had been located at x', then the flux at a position x would be
[Proof:
Hence in general, for S-+S(x1)dx', we find
EXAMPLE: Consider
170 / THE ONESPEED DIFFUSION MODEL OF A NUCLEAR REACTOR THE ONESPEED DIFFUSION THEORY MODEL / 171
for an infinite medium, - oo < x < + oo. Now define G by with boundary conditions such that +(x)< co as x+?co. The solutions to the
homogeneous equation are exp(?x/L), and the particular solution is obviously
--
d2G -LG(~,~~)=--
6 ( x - x') Nx)= S o L 2 / D . Hence we will seek a general solution
dx2 L2 D *
Hence the plane source kernel is just the Green's function for this operator, and we
find
(a) +(;)=o
and a particular solution of the inhomogeneous equation (b) +(-f)=0.
Our boundary conditions require The reader may have already encountered eigenvalue problems in a somewhat
different form:
Adding and subtracting these equations, we find that we must simultaneously where 4, is the eigenfunction corresponding to the eigenvalue A,. However by
require comparing this form with Eq. (5-78), one can easily identify H = d2/dx2, A+ - B 2,
and +tt++tt tx).
Notice that in the example Eq. (5-78) we actually find two types of eigenfunc-
tions: the cosine functions corresponding to odd n and symmetric about the origin,
and the sine functions corresponding to even n and antisymmetric about the origin.
and
Had we restricted ourselves to symmetric sources S(x)= S(-x), we could have
eliminated the antisymmetric solutions [Eq. (5-85)) from further consideration. We
have sketched the first few eigenfunctions for the slab geometry in Figure 5-10.
In acoustics these eigenfunctions would be identified as the normal modes or
How do we achieve this? Certainly we cannot set A , and A, equal to zero since natural harmonics of the system, and this terminology is frequently carried over to
then we would have the "trivial" solution 4(x)=O. Instead we must choose the reactor analysis. Notice that the A, are still undetermined and are, in fact,
parameter B such that these conditions are satisfied. Of course there are many arbitrary. These can be chosen in a number of ways, but for now we will just set
values of B for which this will occur. For example, if we choose A,=O, then any A, = I for convenience.
So now this auxiliary problem has given us an infinite set of solutions cC.,(x).
What good are they? Well, they have a couple of very useful properties. First notice
that the product of any two of these functions will vanish when integrated over the
will give rise to a solution slab unless the functions are identical:
1- 2
dx*.(x)+.(x)=
0, ifm#n
g if m = n .
2'
(5-89)
which obviously satisfies both the differential equation [Eq. (5-78)] and the This property is known as orthogonalify and proves to be of very considerable
boundary conditions. Alternatively, we could have chosen A, =O, in which case usefulness, as we will see in a moment.
The second property of the eigenfunctions +,(x) is that they form a complete set
in the mathematical sense that any reasonably "well-behaved"' function f(x) can
be represented as a linear combination of the +,,(x):
yields solutions
+,(x)=A,sin nnx
n=2,4,6 ,...
7,
a (5-85)
Hence our homogeneous problem can be solved only for certain values of the
parameter B. One refers to the values of B2 for which nontrivial solutions exist to
the homogeneous problem as eigenoalues:
~ , , c o s ( ~ ) n=1,3,5
, ,...
Eigenfunctions: +,,(x) = (5-87)
~ , s i n ( Y ) , n=2,4,6 ,...
FIGURE 5-10. The eigenlunctlonslor a dab lpometrJ.
174 / THE ONESPEED DIFFUSION MODEL OF A NUCLEAR REACTOR THE ONE-SPEED DIFFUSION THEORY MODEL / 175
Of course, such a representation is only of formal interest unless some scheme is Thus we are left with one equation for an infinite number of unknowns, the en.
available for determining the expansion coefficients c,, but this is where orthogon- Fortunately orthogonality once again comes to our rescue. Multiply by J/,(x),
ality~omesin handy. Multiply Eq. (5-90) by $,(x) and integrate over x to find integrate over x, and use orthogonality to find
d -
a -a
2 2 2
Jl"hG)f(d= 2 s /-- , ~ $ m ( ~ ) $ n ~ ) = c m / _ ~ d ~ $ ~ ( x )
2 n-1 2 2
Hence given any function f(x), we can evaluate the appropriate expansion where $,,(x) = cos(nwx/Z), n odd and +,(x) 3 sin( nnx /a'), n even.
coefficients c, by a simple integration. It should be mentioned that for the specific We can rewrite this in a bit more familiar form if we substitute Eq. (5-94) into
example of a finite slab we have been considering, an eigenfunction expansion is Eq. (5-98) and rearrange things a bit:
rimply a Fourier series expansion and is probably already quite familar to most
students.' However, the properties of orthogonality and completeness characteriz-
ing such trigonometric functions also hold for much more general eigenfunc-
tion expansions.
With this background, we are finally ready to return to solve our original
boundary value problem. As we mentioned, the essential idea is to seek the solution
as an expansion in the eigenfunctions qn(x):
Notice that since S(x) is known, we can use orthogonality to determine the source
expansion coefficients [as in Eq. (5-91)] as Note in particular that the Green's function for a finite geometry is no longer a
displacement kernel, that is, a function only of x-x', as it was for infinite
geometries.
This intimate relationship between Green's functions and eigenfunction ex-
~ansionsis actually a very. general result. Suppose we consider the diffusion
~
v Vn+ B"V" ( r )= 0,
Boundary condition: +,,(Fs) -0.
176 / THE ONESPEED DIFFUSION MODEL OF A NUCLEAR REACI'OR
THE ONESPEED DIFFUSION THEORY MODEL / 177
One can demonstrate that these eigenfunctions are orthogonal in the sense that
calculation the heterogeneous nature of the core must be taken into account. One
not only must consider nonuniformities corresponding to fuel pellets, cladding
material, moderator, coolant, control elements, but spatial variations in fuel and
coolant densities due to nonuniform core power densities and temperature distribu-
Since Eq. (5-102) is homogeneous, any solution +,,(r) may be multiplied by a tions as well. Such complexities immediately force one to discard analytical
constant. We will scale each \C;, so that the eigenfunctions are normalized such that methoqs in favor of a direct numerical solution of the diffusion equation. In fact
even when an analytical solution of the diffusion equation is possible, it is
frequently more convenient to bypass this in favor of a numerical solution,
particularly when the analytical solution may involve numerous functions that have
lour earlier slab eigenfunctions +,,(x) would be normalized if we multiplied them to be evaluated numerically in any event, or when parameter studies are required
by ( 2 / ~ ) ' / ~ .The
] set of orthogonal and normalized eigenfunctions (rl,(r)) is said to that may involve a great many such solutions.
be orthonormal. It can also be shown to be complete. Hence we can expand The general procedure is to rewrite the differential diffusion equation in finite
difference form and then solve the resulting system of difference equations on a
digital computer. It is perhaps easiest to illustrate this approach by a very simple
example (sufficiently simple, in fact, to enable analytical solution). Suppose we
wish to solve
where
subject to the boundary conditions characterizing a finite slab of width a :
If we substitute these expansions into Eq. (5-101), we can use orthogonality as (For convenience we will ignore the extrapolation length.)
before to solve for We first discretize the spatial variable x by choosing a set of N + 1 discrete points
equally spaced a distance A = a / N apart (for convenience).
Hence we find
where We now want to rewrite Eq. (5-109) at eath of these discrete points xi, but to do so
we need an approximation for d2+/dx2. Suppose we Taylor expand at xi+, in +
+n(r)+rn (r') terms of its value at the point xi:
G (r, r') = -!-
za I+ L~B,~
as a general result for any geometry (although the eigenfunctions or spatial modes
+,,(r) may be very hard to construct in practice).
dijferenceformula should be a reasonable approximation to the value of d2+/dx2 at subject to boundary or interface conditions that we will leave arbitrary for the
the point xi. moment. Actually we should remark here that this form of the diffusion equation is
If we now use this difference formula to write Eq. (5-109) at any mesh point x,, even a bit too general for most reactor applications. One rarely encounters reactor
we find configurations in which the composition varies in a continuous way from point to
point [i.e., D(x) and Z,(x)]. Rather the system properties are assumed to be
essentially uniform in various subregions of the reactor core (or can be suitably
represented by spatially averaged or "homogenized" properties within each subre-
gion). Hence the far more common situation is one in which the diffusion equation
where again we have defined St=S(x,). We can rearrange this dijference equation [Eq. (5-109)] with constant D, and X,. must be solved in a number of regionsj. We
to rewrite it as will develop the difference equations for the general diffusion equation [Eq.
(5-1 15)], however, since they are not really any more difficult to derive or solve,
and in certain cases they are useful in avoiding technical difficulties arising in less
general approaches (such as the handling of region interfaces).
As in our simple example we begin by setting up our discrete spatial mesh as
shown below, although we will now allow for nonuniform mesh spacing.
Hence we now have reduced Eq. (5-109) to a set of N- 1 algebraic equations for
N + I unknowns (+, $9, 92,...,+N). If we add on the boundary conditions at either There are a variety of schemes that can be used to generate a difference equation
end, say 9, -0, +N-0, we can now imagine solving (or imagine the computer representation of Eq. (5-1 15) on this mesh. We have already considered a simple
solving) this set of algebraic equations. In this particular case, the system of problem in which a Taylor series expansion was used to derive a central difference
algebraic equations can be solved directly using Gaussian elimination." More formula for d2+/dx2. A more common scheme is to integrate the original differen-
generally one must use iterative methods to solve the finite-difference equations. tial equation over an arbitrary mesh interval, and then to suitably approximate
This very simple example illustrates the two essential tasks involved in the these integrals (after an occasional integration by parts) using simple mean values
numerical solution of the diffusion equation: (a) derivation of the corresponding or difference formulas. By way of illi~stration,suppose we integrate Eq. (5-1 15)
difference equations and (b) formulation of a suitable algorithm for solving these over a mesh interval xi - Ai/2 < x < xi + A,+ ,/2 surrounding the mesh point xi.
equations on a digital computer. The methods used will vary from problem to
problem. For example, whereas a direct solution of the difference equations [e.g.,
Eq. (5-1 14)) is possible for one-dimensional problems, iterative methods are re-
quired for two- and three-dimensional problems. Furthermore one generally desires
to work with nonuniform meshes in reactor calculations, to account for the fact
that the neutron flux may vary much more rapidly in certain regions than in others.
In this section we illustrate several of the techniques that are commonly applied
in reactor analysis to the derivation and solution of difference equations. However Let us choose the simplest scheme to approximate the integrals by expressing them
just as in our earlier study of analytical techniques, we do not intend this as the value of the integrand evaluated at the meshpoint xi times the integration
discussion to be a detailed discussion of numerical methods in nuclear reactor interval. For example,
calculations. Instead we refer the interested reader to the extensive literature on
this t ~ ~ i c . ~ " ' ~
The derivative term requires a bit more work. First write depend on x, we return to our earlier results [Eq. (5-1 13)] derived via a Taylor
series expansion. Actually we should note that the coefficients aUdepend only on a
single subscript, j. However double subscripting is useful, for we will rewrite these
algkbraic equations as a matrix equation.
Our final task is to append to these equations two additional equations takirg
into account the boundary conditions. Of course we could simply use the vacuum
To handle d+/dx, we can use a simple two-point difference formula [which can be extrapolated boundary conditions, +,=0, +N = O as before (taking care to place the
derived by subtracting Eqs. (5-1 ION: mesh points xo and xN on these extrapolated boundaries). More general boundary
conditions (such as nonreentrant current) can be developed by taking the final two
difference equations in the set as
and
If we now combine Eqs. (5-1 16), (5-1 17), and (5-120), we arrive at a set of while in spherical coordinates, we find
difference equations very similar to our earlier results
where Hence we can derive difference equations corresponding to these geometries, using
either of the earlier techniques, to find for uniform mesh spacing1'
where now
20
a. . = -
Isr '2
+ X,,
Hence once again we have arrived at a set of A'- 1 three-point difjerence equations
for the N + 1 unknown discretized fluxes, +,+, ...,+N. In the particular case in
which the mesh size Ai is constant and the coefficients D(x) and Z,(x) do not
182 / THE ONESPEED DIFFUSION MODEL OF A NUCLEAR REAmOR THE ONESPEED DIFFUSION THEORY MODEL / 183
where c=0,1,2 for plane, cylindrical, or spherical geometries, respectively. To next task is the determination of a suitable prescription or algorithm for solving this
complete the specification of these difference equations, we will consider the case system of algebraic equations.
of vacuum boundary conditions. For reference we first consider the slab geometry
with uniform mesh spacing A: 3. SOLUTION OF THREE-POINT DIFFERENCE EQUATIONS
Suppose we have developed an appropriate set of difference equations similar
to Eq. (5-121). We must now solve for the discretized fluxes +,.To be more explicit,
let's first write out these equations in detail
,
Notice if we were to apply Eq. (5- 126) to the case i =0 (assuming + - = 0), then this
+,
boundary condition would imply =O. However this is inconsistent with the i = 1
It is useful to rewrite these equations in matrix form
equation. Hence we must ignore the i=O (and i = N ) cases in applying the
difference equation, and use the boundary conditions to eliminate the unknowns +o
and +N from the i - 1 and i = N- 1 equations. For example, the i = 1 equation is
Now consider the case of cylindrical geometry. Our boundary conditions are
now somewhat different. At rN we still have the usual vacuum boundary condition
+N -0. At the origin, we use symmetry to imply
schematically the "forward elimination" on the first two equations in Eq. (5-132):
the matrix into a product of a lower (&) and upper (Y) triangular matrix:I9
As an aside we should observe that while such methods for solving systems of
linear algebraic equations are most easily understood and analyzed (mathemati-
where
cally) in matrix notation, they are most easily programmed when written as a
simple algorithm such as Eq. (5-137). For example, one could simply construct a
loop to generate and store all An and a, using Eq. (5-137) and then evaluate all +,
using Eq. (5- 138).
This algorithm for solving such sets of three-term equations (i.e., inverting
tridiagonal matrices) is easily programmed and executed on a digital computer.
The algorithm would also formally work for solving difference equations
characterizing two- or three-dimensional diffusion problems, however it then en-
We can now substitute back up the matrix to find counters some severe computing limitations. To visualize this more clearly, we will
now briefly comment on the numerical solution of multidimensional diffusion
equations.
4. DERIVATION O F
MULTIDIMENSIONAL DIFFERENCE EQUATIONS
Most detailed neutron diffusion calculations characterizing nuclear reactors
require either two- or three-dimensional treatments. Such details are particularly
important in studying power profiles in large reactors subject to nonuniform fuel
and so on. loading and depletion. Hence we now must consider the numerical solution of the
Thus Gaussian elimination consisting of forward elimination and backward more general d~ffusionequation
substitution can be used to directly solve the difference equations [Eqs. (5-132))
This scheme is particularly important since it frequently appears as an integral part
of the iterative methods used in two- and three-dimensional diffusion problems.
For this reason it is useful to formalize Gaussian elimination a bit by noting that Once again the geometry of interest is discretized into a mesh of cells such as the
what we have in fact accomplished by forward elimination is the factorization of rectangular grids illustrated in Figure 5-1 1. Perhaps the most general way to derive
186 / THE ONE-SPEED DIFFUSION MODEL OF A NUCLEAR REACTOR THE ONESPEED DIFFUSION THEORY MODEL / 187
Our mesh will be defined such that the mesh points are denoted by xo, x,, ..., x i . ..,
FIGURE 5-11. R e d m two- rad Ihree.dlwnslwrl grldc XN; Yo, Y I , ..., y,, ..., y , with mesh spacing Ax and Ay respectively. The most
direct manner in which to derive the difference equations is to use a central
difference formula to approximate
difference equations for the mesh is to integrate the diffusion equation [Eq. (5-143))
over the spatial volume of a given mesh cell, using this to define the spatially
averaged cell properties. In general one can writeM
Here the sum is taken over the adjacent mesh point neighbors j = 1,. .. ,J where
J = 2 , 4, or 6 in I-, 2-, or 3-dimensional Cartesian geometries, while
where the mesh coupling coefficients Iu are determined by the particular mesh
FIGURE $12 l%e twodime~~slonrl
sp.tial mesh
geometry and finite-difference scheme one chooses. For example, in Cartesian
coordinates using essentially the approximation schemes represented by Eqs. (5-
1 16) and (5-120), one would find (We will defer to an exercise at the end of the chapter the more general derivation
for nonuniform media using mesh cell integration.) Substituting these expressions
into Eq. (5-152) evaluated at the mesh point (xi,yj)we find
where we define
The difference equations representing Eq. (5-143) then take the form
where i runs over all of the mesh points. As before, we can use boundary conditions to specify + o J , + N j , + i , o , and $J~,,.
1811 / THE ONESPEED DIFFUSION MODEL OF A NUCLEAR REA(;TOR THE ONESPEED DIFFUSION THEORY MODEL / 189
5. ITERATIVE SOLUTION OF sweep reduced the original tridiagonal matrix to a form with only two diagonals,
MULTIDIMENSIONAL DIFFERENCE EQUATIONS while the backward substitution step completed the matrix inversion
We now turn our attention to the solution of these difference equations. Our
first task is to cast the set of equations into matrix form. This requires first
assigning a single index to each mesh point. For example, in a two-dimensional
mesh we could label the mesh points as
For a two-dimensional problem, the matrix structure takes the following form-in After a bit of examination, it becomes apparent that when a similar forward sweep
this case, for a five-by-four mesh point array is conducted on the five-diagonal matrix characterizing two-dimensional problems,
the result is to fill in all of the zero entries between the main and outer diagonal.
This implies that one will require considerably more computer memory to allow
such a direct inversion of the matrix. Such a direct algorithm is also rather
complicated to program and leads to problems resulting from computer round-off
I
\ \-(-3M-"IH.7+% +* + ~ -3-
+
-32. ~ error. For these reasons, it is far more efficient to use an iterative procedure to
invert such matrices when N is large, since such schemes attempt to preserve the
( a l I) (AT 01 (A,,, A , u7,d2 (A,. ,,12X (A7..12H r ,. ) 2 k 7 ) =(.I)
sparse structure of the original matrix in their operations.
(5-155) Let's illustrate the basic idea with a simple example: Suppose we wish to invert a
Notice that the tridiagonal form we encountered in the one-dimensional case has matrix &--that is, we wish to solve
now been augmented by two additional side-band diagonals-as we might have
expected, since on closer examination we find the two-dimensional case yields a
five-point difference equation. We first compose A_ into its diagonal and off-diagonal elements
Similarly, assigning a single index to each mesh point of a three-dimensional
problem yields the matrix structure which corresponds to a seven-point difference
equation.
Now let's consider how we might solve such systems-that is, invert such
matrices. Since Gaussian elimination can be applied to any matrix (formally, at
least), we might first consider applying this technique to obtain a direct inversion.
Recall that for a one-dimensional diffusion equation, the forward elimination
190 / THE ONESPEED DIFFUSION MODEL OF A NUCLEAR REACTOR THE ONESPEED DIFFUSION THEORY MODEL / 191
Hence suppose we use Eq. (5-160) to rewrite Eq. (5-159) first as explicitly in terms of the algebraic system
ep=,Bp+s
+ a12+Jm)+ aI3+$")+. . . + alN+Am) = S l
and then invert 4 to find
a2,+jm'+ . . . + aZN+km) = S2
Now is where the iterative philosophy comes in. Suppose we guess 2 on the
right-hand side--call the guess (O)-and then use it to calculate a new guess, aNI+lm) + aN2+!jm)
+ + - .. +
aN3+Jm)
as
(5-167)
Hence we can solve for the m + I flux iterate immediately as
It should be noted here that the Jacobi scheme does not use all of the available
Hopefully, then, as m becomes large, we converge to the true solution information during each iteration. For example, if the equations are solved in
sequence from i = I to i = N, as they would be on a computer, then the solution of
the first equation yields +$'"+I); but to find +jm+') using the second equation, +im)
is used rather than the improved estimate +Im+'). Similarly, solving the third
) +4"+') which
equation for +jm+') makes use of +fm)and +im)rather than + [ ' " + Iand
are known. If these estimates are used as soon as they are generated, a more
efficient iterative scheme known as the Gauss-Seidel or successive relaxation method
Hence the general idea behind such iterative schemes is to generate improved is obtained. In this case, the system of equations in each iteration is solved as
guesses or iterates ( " ) by solving the original system of equations in an approxi-
mate, but efficient, manner. We continue such an iterative process until two
successive iterates +('") and + ( " + I ) are sufficiently close together, at which point
the iteration is stopped a n h *("'+I) is regarded as the solution. Notice that
throughout the iterative process, we maintain the sparse structure of the original
five-diagonal matrix 4,thereby significantly reducing storage and calculational
requirements.
The particular scheme we have presented is known as the Jacobi-Richarchon or
PointJacobi method, and although it is a very simple scheme, it has the drawback
that it converges very slowly. The reader might very roughly think of the conver-
gence rate of such iterative processes as being determined by how big a chunk of
the original matrix he is willing to invert on each iteration. (More precisely, the
convergence rate is determined by the size of the matrix norm of D - l L . ) In the
Point-Jacobi method, only a relatively small bit of the matrix, its main diagonal, is and the solution is
inverted on each step and hence we might expect convergence to be slow. (As an
extreme example, one bites off a much bigger chunk in Gaussian elimination-the
whole matrix A_--and hence only a single iteration is needed.) One can accelerate
this convergence in several ways. First, one could attempt to invert a bigger chunk
of A on each iteration. It is also possible to use information about the next flux
.
From solution of
/
From previous (5-170)
iterate during an iterative step. Finally, one can extrapolate from earlier flux previous equations (mth) iteration
it'erates in order to more rapidly approach the true solution. in current ( m + I)
To understand how to improve the Jacobi iterative scheme, let's write i t out iteration
192 / THE ONESPEED DIFFUSION MODEL OF A NUCLEAR REACTOR THE ONESPEED DIFFUSION THEORY MODEL / 193
This can be rewritten in matrix form by decomposing 4 into the sum of an upper
and lower triangular matrix:
k-
SOR @2w I..C
dL
SR @jm'
Previous ,,I)
Here L contains elements of the main diagonal and below it, while Y contains iterate @i
N
- au+jm)]
t
Gauss-Seidel From SOR in
'-i+i/
From SOR in
(5-174)
estimate current iteration previous iteration
[(m+ I)stl [mthl
Now +/'"+') is calculated as a linear combination of +i ("+i ) and the previous SOR
iterate Again iterative methods are utilized in which the outer diagonal elements are
handled in a manner similar to those used in two-dimensional problems. However
there is some reduction in iterative convergence rates due to a loss of procedure
implicitness caused by the additional diagonal elements.
Here the extrapolation or acceleration parameter w ranges between I and 2. Of Such iterative algorithms for the solution of the finite difference equations
course for w = I we return to the Gauss-Seidel method in which no extrapolation is characterizing two- or three-dimensional diffusion problems are frequently referred
used. The iterative algorithm for each element can then be written as to as inner iterufions. This terminology arises from the fact that in nucle~rreactor
criticality calculations, the solution of the diffusion equation is itself imbedded in
yet another iterative scheme-the so-called outer or source iterations-necessary to
handle the presence of a fission term. We will study this latter scheme in Section
5-IV.
194 / THE ONESPEED DIFFUSION MODEL OF A NUCLEAR REAfTOR THE ONESPEED DIFFUSION THEORY MODEL / 1%
6. NODAL METHODS 5-14. We now integrate Eq. (5-178) over the volume Vnof the nth node cell
There are many instances in nuclear reactor analysis in which one requires a
full three-dimensional calculation of the neutron flux, for example, in core fuel
depletion or control rod ejection studies. Although a direct numerical solution of
the diffusion equation can be performed on a modern digital computer, it is
extremely expensive to do so, particularly when a series of such calculations would
be required for a parameter study. We desire a scheme for determining the
three-dimensional core flux distribution that avoids the large storage and execution
time requirements of a direct finite difference treatment of the diffusion eauation. If we define the spatial averages over the nodal cells,
a
Such scheme is provided by so-called nodal r n e t h o d ~ . ~ "The
~ ~ general idea is to
decompose the reactor core into relatively large subregions or node cells in which
the material composition and flux are assumed uniform (or at least treated in an
average sense). One then attempts to determine the coupling coefficients
characterizing node cell to node cell leakage and then to determine the node cell
fluxes themselves.
To develop this approach in more detail, consider the neutron diffusion equation
in its general time-independent form given by Eq. (5-143). Now we know that we
can formally write the solution to this equation as
and allowing neutron transfer to only the six nearest neighbor cells. The transfer
coefficients are represented as linear combinations of several simple trial functions
(e.g., from one-dimensional slab geometry calculations). The blending coefficients
in this representation are then determined by comparison with more accurate finite
Scatter
difference benchmark calculations and lots of experience, fiddling, and fudging.
If the transfer coefficients are properly chosen, then such nodal methods can be I I
extremely useful in generating three-dimensional flux distributions when only
limited accuracy is required. Unfortunately such empirical schemes for choosing
the K,,. are quite problem-sensitive and require a good deal of experience on the
part of the reactor analyst.
We will leave numerical methods for solving diffusion equations until later when
we must generalize these methods to account for fission processes and energy-
dependence. The above discussion has been an admittedly curse description of I absorption
numerical methods for solving differential equations. There is a vast literature on
this subject that provides the details of the methods we have so briefly ~utlined.'~-'~
And perhaps the most valid argument for presenting only a brief sketch of such
numerical methods lies in the recognition that these topics are of such vital
importance to the practicing nuclear engineer, that he almost certainly will have
had or will take further courses on numerical analysis in any event.
Scatter - Leakage @/
fission^
II -
Resonance
Fission J
--
a x = 0 -
a
2 2 FIGURE 5-17. 'Ibe slab re&m
FIGURE 5-16
neutron flux in such a reactor is
Neutron prortllcr h o l v e d la a f a t mcta.
If this is the only source of neutrons in the reactor,+ then the appropriate diffusion
equation becomes
-(X ) (symmetric),
with initial condition: +(x,O) = +,(x) = $ J ~
If we substitute this form into Eq. (5-189) and divide by #(x)T(t), we find
C. The Time-Dependent "Slab" Reactor
I . GENERAL SOLUTION
We will begin our study of nuclear reactor behavior as described by the
one-speed diffusion equation by considering a uniform slab of fissile material Here we have noted that since we have a function only of x set equal to a function
characterized by cross sections Z,. 8,, and Z,. This unrealistic appearing "slab only of r, both terms must in fact be equal to a constant. We have named this
reactor" is chosen to introduce many of the concepts of nuclear reactor analysis, constant - A . However A is as yet unknown.
since its one dimensional geometry greatly facilitates the detailed solution of the Hence the separation of variables given by Eq. (5-190) has reduced the original
one speed diffusion equation.* The appropriate mathematical description of the partial differential equation in two variables to two ordinary differential equations:
t Actually we should hedge here a bit. Equation (5-187) actually represents only "prompt"
neutrons, that is, those born instantaneously in fission. The "delayed" neutrons arising from
fission product decay require a slightly different treatment. We will defer this modification until
Chapter 6.
In this sense it is somewhat akin to the "vibrating string" or "simple harmonic oscillator" We can easily solve the time-dependent equation
problems in physics that also get beaten to death since they contain most of the interesting
physics-and yet are easy to solve.
ZW / THE ONESPEED DIFFUSION MODEL OF A NUCLEAR REACTOR
THE ONESPEED DIFFUSION THEORY MODEL / 201
where T(0) is an initial value which must be determined later. To solve the
space-dependent equation, we must tack on the boundary conditions: factor:
d
where we have noted that the expansion coefficients now must be time-dependent,
If we identify Eq. (5-194) as the same problem, it is apparent that we must choose then we could have immediately arrived at an equation for the T,
These values of A, are known as the time eigenwlues of the equation, since they
characterize the time decay in Eq. (5-193). The general solution to Eq. (5-189) must by substituting Eq. (5-203) into the original equation (5-189) and using the
therefore be of the form orthogonality property of the eigenfunctions. This alternative approach is
frequently useful when encountering problems in which sources are present, be-
cause the separation of variables approach we have presented applies only to
+(x,t)= A, exp(-X,t)cos ?. (5- 198) homogeneous equations.
n
odd Finally, note that although we initially sought solutions +(x)T(t) which were
separable in x and t, these solutions were eventually superimposed to yield a
This solution automatically satisfies the boundary conditions. To determine the A,, nonseparable function of space and time [cf. Eq. (5-201)]. Hence separation of
we use the initial condition to write variables does certainly not imply a separable solution. Interestingly enough,
Notice that we have qualified this definition by specifically demanding that the flux
be time-independent in the absence of a source. As we have seen in Chapter 3 (and
will see later in the problem set at the end of this chapter), a source present in a
critical system will give rise to an increase in flux that is linear in time.
If we write out the general solution for the flux
m
+(x,t)=A,exp(-A,t)cosB,x+ Anexp(-A,r)cosB,x, (5-209)
n-3
n odd
it is evident that the requirement for a time-independent flux is just that the
fundamental time eigenvalue vanish
fundamental mode shape. Of course, we have implicitly assumed that A, will not be since then the higher modes will have negative A, and decay out in time, leaving
zero. The coefficient of the fundamental mode is just just
+(x, r)-+A,cos B,x # function of time. (5-21 1)
3. THE CRITICALITY CONDITION In particular notice that by increasing the core size we decrease B ,: while by
increasing the concentration of fissile material we increase 1,and hence B:. Both
Let us now see what is required to make the flux distribution in the reactor of these modifications would therefore tend to enhance core multiplication.
time-independent-that is, to make the fission chain reaction steady-state. We will Yet recall that in Chapter 3 we expressed the criticality condition in terms of the
define this situation to be that of reactor criricaliry: multiplication factor k. We can make the connection between these two criteria if
Criticality r when a time-independent neutron flux can be sus- we write the time eigenvalue as
tained in the reactor (in the absence of sources
other than fission).
2Q4 / THE ONESPEED DIFFWSION MODEL OF A NUCLEAR REACTOR THE ONESPEED DIFFUSION THEORY MODEL / UH
Now recall that (ox)-' is the mean lifetime for a given neutron-nuclear reaction to then it is apparent that the various forms of the criticality condition are indeed
occur. Hence (ox,)-' must be just the mean lifetime of a neutron to absorption equivalent:
(ignoring leakage-that is, in an infinite medium). Furthermore, we can identify
In particular notice that by using PNL=(l+ L28,2)-I we can avoid the analysis of
the initial value problem and proceed directly to the criticality condition
Now the only remaining task is to identify (1 + L2B;)-I. Recall that the rate of
neutron leakage is given by
rate L
Leakage = dS .J = d + v . J = - dlr ~ v 2 + ,
Jv Jv
(5-2 18)
We will return later to consider how these results can be applied to reactor
criticality studies, but first we will extend them to more general reactor geometries.
where we have used both Gauss's theorem and the diffusion approximation. Hence
we can write D. The Criticality Condition for More General Bare Geometries
Note that the only quantity characteristic of the reactor size or geometry that
appears in k or PNLis the geometric buckling, 8;. For the case of a slab reactor of
Rate of neutron absorption = J$3r~a+ width a we found ~ l = ( n / a ' ) ~We
. might suspect that for more general geometries
Rate of neutron absorption plus leakage
Jv d3rXa+- 1"d3rDV2+ we need only replace this by the geometric buckling characterizing the specific
geometry under consideration. This suspicion is in fact easily verified, but only for
so-called "bare" geometries in which the reactor composition is uniform. For the
more complicated multiregion geometries, such as reactors composed of a core
surrounded by a reflecting material, one can no longer derive simple expressions
for PNLor k in terms of the reactor geometry and composition.
Consider, then, a bare reactor of uniform composition surrounded by a free
However we can identify this ratio as just nonreentrant surface characterized by vacuum boundary conditions. If the reactor
is critical then the neutron flux must satisfy the steady-state diffusion equation
Nonleakage probability EE PNL=- 1 (5-220)
I + L'B,~ '
Therefore we can interpret subject to the boundary condition +(i,)=O for is on the extrapolated surface. Of
course in general there will be no solution to this equation unless we have
1 -,=neutron lifetime
(5-22 1) happened to hit on just the right combination of composition and system size.
( v ~ . ) ( I+ L~B:) 'P ~ ~ ( & ) = in a finite reactor, To see this more clearly, divide Eq. (5-227) by - D so that it can be written as
since we have just reduced the lifetime to absorption in an infinite medium to take
account of neutron leakage. If we now combine Eqs. (5-217) and (5-220), we find
that the multiplication factor k for this model becomes just
boundary condition: +(is) = 0.
Sometimes Eq. (5-228) is written in a somewhat different form as
Thus we can identify our fundamental time eigenvalue as just the inverse of the
reactor period boundary condition: +(is) =O.
-A - k-1 1 Now notice that this equation is identical to that which generates the spatial
1 - 7 ' 7
eigenfunctions for this geometry
If we also recall from Eq. (5-202) that
VF, + B,F.,(r) = 0, (5-230)
boundary condition: +(is) = 0.
206 / THE ONE-SPEED DIFFUSION MODEL OF A NUCLEAR REACTOR THE ONESPEED DIFFUSION THEORY MODEL / 207
We know that this latter equation has nontrivial solutions h ( r ) only for certain This is just the local thermal power density at position r in the core. Hence the total
values of the parameter B 2, the eigenvalues 8.: Hence by comparing Eqs. (5-228) power generated by the core is just the integral of the power density over the core
and (5-230) we find that the steady-state diffusion equation for the flux Hr) will volume
only have nontrivial solutions when the core composition is such that (vX, - Z,)/D
is equal to an eigenvalue B,: and then the flux Nr) will be given by the
corresponding eigenfunction +,,(r).
However since there are an infinite number of possible eigenvalues Bi, we might This relation can be used to determine the magnitude of the flux in terms of the
be tempted to think that there are an infinite number of values (~2,-Z,)/D r B: core thermal power level.
for which the reactor is critical. However it should be recalled that in the case of a Thus we now have developed a rather simple scheme to study the criticality of a
slab geometry, only the lowest eigenvalue B:=(T/~')~GB,' had a corresponding nuclear reactor-at least a bare, uniform reactor. The only mathematical effort
eigenfunction $,(x) -cosnx/d that was everywhere positive. The eigenfunctions or required is the solution of the Helmhotz equation characterizing the geometry of
spatial modes & ( x ) correspbnding to higher ;igenvaiues oscillated about zero. This
interest for the geometric buckling 8: (the fundamental eigenvalue B:) and the
same feature also characterizes the eigenfunctions +,,(r) of more general geometries.
critical flux shape +(r) [the fundamental spatial eigenfunction +,(r)]. To illustrate
Only the eigenfunction +,(r) corresponding to the smallest eigenvalue B: is
everywhere nonnegative. Since the neutron flux can never be negative, it is how these quantities are determined, we will consider a simple yet very important
apparent that the only solution to the eigenvalue problem Eq. (5-230) physically example:
relevant is that corresponding to the smallest eigenvalue, B:= 8.: Hence for the
EXAMPLE: A Right Circular Cylindrical Core
reactor to be critical we require
The most common reactor core shape is that of a right circular cylinder of height
H a n d radius R. (Actually a sphere would be the more optimum geometry from the
aspect of minimizing neutron leakage, but spheres are very inconvenient geometries
to pass coolants through.) The appropriate form of the Helmholtz equation is then
just as for the slab.
Thus we can continue to use P,,=(l+ L'B,')-I and k=(vZ,/Z&(I + L2B,')-I
for more general bare geometries provided we identify the geometric buckling B
:
as the smallest eigenvalue B: of the Helmhotz equation subject to boundary conditions
subject to the boundary conditions that +(r) vanish on the extrapolated boundary
of the reactor. The corresponding critical flux distribution +(r) is then given by the Since this is a homogeneous partial differential equation, we can seek its solution
fundamental eigenfunction +,(r), which is everywhere nonnegative. using separation of variables
It should be pointed out that although the Helmhotz equation (5-232) will
provide us with the flux shape in a critical reactor, it will tell us nothing about the
magnitude of the flux. Since it is a homogeneous equation, if +(r) is a solution, then Then if we substitute this form into Eq. (5-235), we arrive at two ordinary
any multiple of +(r) is also a solution. Of course the magnitude of the flux was
determined for us by the initial condition +o(r) when we studied reactor criticality
by solving the full time-dependent diffusion equation (5-187). However this latter
approach is far too cumbersome to use in practice.
Instead we merely note that a critical reactor can operate at any flux level-at
least mathematically. (Of course one must provide for adequate core cooling,
shielding, etc., but these factors are extraneous to our present model of the reactor
so we won't worry about them here.) Hence we will merely assume that the
magnitude of the flux is determined by the desired thermal power output of the
core. If the usable energy produced per fission event is w,, then the thermal energy
deposited in the core per unit volume per second is just given in terms of the fission
reaction rate density as
FIGURE 3-19. Flolte cylindrical reactor core
208 / THE ONESPEED DIFFUSION MODEL OF A NUCLEAR REACTOR THE ONESPEED DIFFUSION THEORY MODEL / 209
differential equations
Thus we find
It should perhaps be mentioned that since reactor cores are fabricated from either
square- or hexagonally-shaped fuel assemblies, one can only approximate such
where the separation constants 'a and A2 are constrained by the relationship cylindrical geometries. However for most purposes one can assume the reactor core
+
B'- a 2 A'. Each of these equations represents a separate eigenvalue problem that is essentially a right circular cylinder.
one can use to determine u and A (and hence 8'). The eigenfunctions and One can proceed in a very similar manner to analyze other bare core geometries.
eigenvalues of the axial equation are well known to us For convenience, we have tabulated the geometric buckling and critical flux profile
in other common geometries in Table 5-1.
Slab
Sin? Yo(ar)+m as r-0, we must set C r O . Applying our boundary condition at
r = R, we find
Infinite
Cylinder
where v, are the zeros of J,. In particular, the smallest such zero is ~ ~ ~ 2 ....
405
(kind of like n to a Bessel function). Hence we find the eigenfunctions and
eigenvalues generated by the radial equation (5-236a) are just
Sphere (2)'
Therefore, consistent with our prescription of seeking the smallest value of B2 as
our geometric buckling, we find
Finite
corresponding to a spatial flux shape Cylinder J0 -) cos( f )
will yield a critical reactor. For example, if the material composition is specified, It should be noted that the transport cross sections used in this example have
one can compute the material buckling B: in terms of the macroscopic cross been artifically adjusted (reduced by almost an order of magnitude) to take some
sections using Eq. (5-213). Then by using the criticality condition B= : B: (along account of fast neutron leakage which would normally not be described by a
with Table 5-I), one can infer the core dimensions that will yield a critical system. one-speed model.
In the more usual situation the nuclear designer will be given B i (rather than B ): We can use these cross sections to calculate a number of important parameters
since the core dimensions are determined by limitations on core thermal perfor- characterizing the PWR core:
mance (and not nuclear considerations). That is, the core must be built--a
sufficiently large size to avoid excessively high temperatures for a desired power Diffusion coefficient: D -9.21 cm
output. The nuclear designer must then determine the fuel concentration or loading Infinite multiplication constant: k, = vZ,/Z,= 1.025
(i.e., B;) that not only will result in a critical system, but will also allow the core to
operate at a rated power for a given time period. Material buckling: B;=(VZ,-Z,)/ D =4.13 x cm-2
EXAMPLE: As a specific example we will study the one-speed diffusion model Extrapolation distance: zo= 0.7 I&= 19.6cm
of a bare, homogeneous cylindrical reactor with material composition representa- Leakage fraction for a critical core: 1 - P,, -0.025 (5-245)
tive of that of a modern PWR such as described in Appendix H. We will use
"homogenized" number densities corresponding to a PWR core operating at full
Next we will compute the critical core dimensions. If we assume that the core
power conditions and containing a concentration of 2210 ppm of natural boron (as
height is fixed at 370 cm by thermal considerations, then we can determine the
boric acid) dissolved in the water coolant for control purposes. The fuel is taken as
radius at which a core with such a composition will be critical. First calculate the
UO, enriched to 2.78% ,"u. In Table 5-2 we have listed the number densities and
axial buckling 8:
microscopic one-speed cross sections for this core. (Here the cross sections are
actually averages over the neutron energy distribution in such a reactor core.)
This data can be used to calculate the macroscopic cross sections tabulated in
Table 5-3. Here we have also included the relative absorption rates in each material (It should be noted that this is somewhat smaller than the radius of 180 cm for a
which serve as a measure of neutron balance within the core. typical PWR core. This illustrates the limitation of such a one-speed model for
obtaining quantitative estimates in reactor analysis.)
TABLE 5 3 M.cmrcaple C r m Secths
while in the reflector we will seek a solution satisfying the vacuum boundary
condition (c)
where the reflector diffusion length LR.=(DR/Z,R)+.We now apply the interface
boundary conditions to find
d2C#sC
Core: -DC-+p:-vx:)C#sC(x)-0, 0<X<f,
dx Notice that this equation represents a relation between reactor composition
Reflectoo -D "*
+ -2 R
dx2
z ~ c # ~ ~ (0,
x) Za < x < ' + i ,
2
(5-249) (D C, B:, D R,LR) and size ( a ,b) that must be satisfied if a solution to the steady-
state diffusion equations (5-249) is to exist. Hence this is just the reactor criticality
condition for this particular geometry. Admittedly, it doesn't look anything like our
subject to the set of boundary conditions: earlier condition, B:= B:, that characterized a bare reactor. In fact the criticality
condition for a reflected reactor is transcendental-one cannot obtain an explicit
solution for the critical size or composition. Instead, either numerical or graphical
techniques must be used. The latter technique is more useful for our present
discussion. Rewrite Eq. (5-254) as
Note that we have used the reactor symmetry to narrow our attention to the range If we plot the LHS against (~:a/2), we can then determine the solution of this
of positive x . As we noted earlier in Section 5-111-D this problem will have no transcendental equation graphically by noting where it intersects the value of the
solution unless we choose the proper combination of core composition and size. RHS, as shown in Figure 5-21. [Actually since there will be many such intersec-
. We would anticipate that a criticality condition relating these core characteristics tions, we are only interested in the lowest value of (B;a/2).] From this graph we
would emerge in the course of our analysis. notice that the critical value of B: must be such that
The general approach, as always, is to defermine the general solutions in the core
and reflector and then use the boundary conditions to determine the unknown
coefficients. In the core the general solution will be
in contrast to the bare (unreflected) core in which B,c2=(n/q2. Hence we see that
the width o required for criticality is somewhat smaller when a reflector is added,
where we have utilized the symmetry of the core to discard the sine term. Here the but we would expect this since a reflector is added primarily to reduce neutron
material buckling characterizing the core is defined by leakage.
It is conventional to define the difference between bare and reflected core
dimensions as the reflector savings 6:
S = a (bare) - a (reflected). (5-257)
114 / THE ONESPEED DIFFUSION MODEL OF A NUCLEAR REACXOR THE ONESPEED DIFFUSION THEORY MODEL / 215
For example, the reflector savings for our slab reactor can be written as
Then for some value of k, we assert that this equation will always have a solution.
The idea is to pick a core size and composition and solve the above equation while
determining k. If k should happen to be unity, we have chosen the critical size and
composition. If k # I, however, we must choose a new size and composition and
For a thick reflector bW LR this simplifies to repeat the calculation. As one might expect, k turns out to indeed be the multipli-
cation factor we defined earlier in Chapter 3, as we will demonstrate later.
We could give a formal mathematical proof that Eq. (5-262), or its generaliza-
tions will always have a solution for some k, but it is more convenient to simply
argue physically that since varying k will vary the effective fuel concentration
which is essentially a measure of the maximum reflector savings that can be NF+NF/k, one can always achieve a critical system by making k sufficiently
realized. small.
Reflectors serve another function besides reducing neutron leakage. They tend to Sometimes a slightly different formulation is used in which one pretends that v,
flatten the flux and hence the power distribution in the reactor core. Unfortunately the number of neutrons emitted per fission, is in fact variable. (Of course it isn't,
one-speed diffusion theory is not adequate to describe this effect, which results in a but it is a useful device to regard it as adjustable for the moment.) Now physically
peaking of the thermal flux in the reflector region (see Figure 5-19), so we must we know that there must be some value of v , call it v, that will yield a nontrivial
defer a further discussion of reflected cores until we have developed multigroup solution to
diffusion theory.
IV. REACTOR CRITICALITY CALCULATIONS regardless of what composition or geometry we have chosen. Hence the idea is to
determine this vc, then readjust composition and geometry until we have forced
A. Introduction
Let us now turn to the very important topic of determining the composition
or size of a reactor that will yield criticality. It should be apparent after the last
example that in most practical reactor designs one cannot simply determine the If we compare this approach to our earlier scheme in which we calculate k, it is
geometric buckling for a core geometry and then use B: = B: to arrive at critical- evident that
ity.
One "brute force" procedure would be to determine the lowest time eigenvalue
of
From a mathematical point of view, each of these approaches introduces a new
parameter into the steady-state diffusion equation, either k or v, which can then
be regarded as an eigenoalue in a subsequent analysis. Once this eigenvalue has
been calculated, one can return and readjust composition and geometry in an effort
and then keep adjusting things until A-0. However this is rather-awkward, and to force this eigenvalue to a desired value (e.g., k + l or v,+v). Hence the criticality
216 / THE ONESPEED DIFFUSION MODEL OF A NUCLEAR REACTOR THE ONE-SPEED DIFNSlON THEORY MODEL / 217
calculation is converted into a sequence of eigenvalue problems for the criticalily using our earlier procedures. With this solution, we can now explicitly calculate the
eigenwlue (sometimes also called the multiplication eigenoalue). fission source resulting from this flux + ( I ) as
Of course, in general there will be a set of criticality eigenvalues k, correspond-
ing to the eigenvalue problem represented by Eq. (5-262). For example, our earlier
analysis for the bare slab reactor indicated the existence of the set of eigenvalues:
k,-(vX,/Zb(I + L28,2)-' where ~ : = ( n r / Z ) ~ ,n = 1,2 .... Only the largest such This can then be taken as a new estimate of the fission source and used to generate
eigenvalue (in this case, k , ) will correspond to an everywhere-nonnegative flux a new flux, +(2), and so on-provided we can also generate improved estimates of
distribution Hr) and hence to a critical reactor configuration. We will refer to the k. That is, we can iteratively solve for an improved source estimate S ( " + ' ) from an
largest criticality eigenvalue k , as the effectiw multiplication factor and denote it by earlier estimate S(") by solving
kep As we will see below, ken can be identified as the multiplication factor k for the
reactor core defined earlier in terms of fission neutron generations in Chapter 3.
We can now use this relationship to compute a new guess of k ( " + ' )from + ( " + I ) and
k(").
218 / THE ONESPEED DIFFUSION MODEL OF A NUCLEAR REACTOR THE ONESPEED DIFFUSION THEORY MODEL / 219
We should note that this prescription is quite consistent with our earlier in-
terpretation of k as the multiplication factor-that is, the ratio of the number of
neutrons in two consecutive fission generations-if we note that a factor of kc")
must be inserted in the denominator since [k(")]-'S(") is in fact the effective fission
I Guess core geometry
and composition
-
Criticality
search
Notice that by scaling the source term appearing in the diffusion equation (5-266)
by a factor of Ilk(") in each iteration, we will prevent the rapid growth or decrease
of successive source iterates (causing possible overflow or underflow) in the event iterations
that a number of iterations are required when k is not close to unity. That is,
dividing the source term by k(") removes the dependence of the flux iterate $("+I)
on n [at least as +("+I) approaches the true solution].
This iterative scheme to determine the effective multiplication factor kd and the
corresponding flux Hr) is known as the power iteration or source iteration method.
The iterations themselves are known as outer or source iterations.
In addition to such outer or source iterations, one will also be required to
perform inner iterations to solve the diffusion problem
C. Source Extrapolation
Needless to say, there is strong incentive to perform as few iterations as
possible in converging to the desired accuracy. For that reason, one usually
attempts to accelerate the source iteration convergence by extrapolating ahead to a
new source guess. This is accomplished by introducing an extraiolati&n parameter
FIGURE 522.
2, FINISHED!
Now suppose we were to uniformly modify or perturb the absorption cross section ing a function by vXf(r).If we write
to a new value
( f , ~ g ) = d 3 r ~ v Z f g J=: 3 r ( v ~ f f ) * g
where we will assume that the perturbation SX, is smalCthat is, = ( v X , f , g )= ( ~ % g ) , (5-286)
where we have merely shuffled Zf(r) around in the integral (noting that 2, is real)
to identify
Then the value of k' corresponding to the perturbed core can be written as
Notice that in this case, F t and F are in fact identical. We refer to such operators
as being self-ac$oint.
where we have expanded k' in SZ,/Z, and have neglected all terms of higher than For a more complicated example, consider the spatial derivatives in the diffusion
first order in the perturbation ( 8 x 3 . This has allowed us to express the perturbed operator M :
multiplication factor k' in terms of the unperturbed multiplication k and the
perturbation 6 2 , .
These general features appear in applications of perturbation theory to more
general problems in nuclear reactor analysis in which the perturbations may be
localized or in which the multigroup diffusion equations are used as the basic Now if we use the vector identity
model of the core behavior. Although the general ideas are essentially as simple as
those in the example above, it is necessary to introduce a few mathematical V.ab=aV.b+b.Va,
preliminaries. (For more details, the reader is referred to Appendix E.)
We will describe the multiplication of the core by the criticality eigenvalue we can rewrite this as
problem [Eq. (5-266))
( f , ~ . D v g ) v= ( d ~ r ~ . [ f * D( vd ~3 ]r -[ ~ r . D v g ] . (5-290)
Using Gauss's law, we can convert the first term into an integral over the surface:
where we will leave it as understood that the solution of this equation, +(r), must
satisfy appropriate boundary conditions such as +(F,)=O on the surface of the core.
Now suppose we define the inner product ( f , g ) between any two functions f ( r ) However since we require that j and g vanish on the surface, this term vanishes. If
and g(r) as we repeat this procedure we find we can rewrite
where f*(r) denotes the complex conjugate off (r), and Y is the core volume.
We can now use this inner product to define the operator ~t adjoint to the Hence we find that
operator M as that operator M t for which
V~D~O=V.DVO
define the adjoint flux +t as the corresponding solution of we can use Eq. (5-300) to find the change in reactivity due to the addition of an
absorption cross section SZa(r)
Note here that the perturbation in the core absorption appears as a perturbation
6M in the diffusion operator Then neglecting second and higher order quantities in the perturbation-that is,
using first order perturbation theory-we find
To calculate the change in k , first take the scalar product of Eq. (5-296) with the
adjoint flux +t characterizing the unperturbed core, that is, satisfying Eq. (5-294),
Since the one-speed diffusion operator is self-adjoint, we know $t=+, and hence
we find
Now using the definition Eq. (5-285) of the adjoint operator, we find
Hence we find
It should be noted that all of this analysis was exact until we neglected second-
order terms in Eq. (5-304). Thus, we have calculated a first-order estimate of the
reactivity change Ap due to introducing a localized absorber 8Z,(r) in terms of the
unperturbed flux distribution.
We could now calculate 6k = k'- k. However it is far more convenient to define
the core reactiviry EXAMPLE: Consider a bare slab reactor characterized by one-group constants
D, Z, and vZ,. We will perturb this reactor by imagining that an additional
absorber is uniformly inserted in the region O < x < h. One might consider this to be
a model of a bank of control rods inserted to a depth h in the core. Of course to
allow the application of perturbation theory, we must assume this absorption to be
which essentially measures the deviation of the core multiplication from unity. relatively small.
Then since the perturbation in reactivity is just Hence our perturbation is
224 / THE ONESPEED DIFFUSION MODEL OF A NUCLEAR REAtXOR THE ONESPEED DIFFUSION THEORY MODEL / 225
Consider the adjoint problem 7. K. M. Case and P. F. Zweifel, Linear Transport Theory, Addison-Wesley, Reading,
Mass. (1967).
8. G. Arfken, Mathematical Methods for Physicists, 2nd Edition, Academic, New York
(1970), pp. 733-739.
(Of course for one-speed diffusion theory, M = M t is self-adjoint, but we will retain 9. Ibid., pp. 748-768.
10. E. A. Coddington, Ordinary Differential Equations, Prentice-Hall, Englewood Cliffs, N.
the generality for a bit.) Notice that we have allowed the source St(r) appearing in
J. (1961). p. 67.
the adjoint equation t o differ from that in Eq. (5-316). 11. B. Carnahan, H. A. Luther, and J. 0.Wilkes, Applied Numerical Methoak, Wiley, New
Now suppose we multiply Eq. (5-316) by gt a n d integrate over r, then multiply York (1969).
+
Eq. (5-317) b y a n d integrate, a n d then subtract these two results to find 12. M. Clark, Jr., and K. F. Hansen, Numerical Methods of Reactor Analysis, Academic,
New York (1964).
13. E. Wachspress, Iteraticv Solution of Elliptic Systems, Prentice-Hall, Englewood Cliffs, N.
J. (1966).
However by the definition of the adjoint Mt- Ft, the LHS is zero. Hence we find: 14. R. S. Varga, Matrix Iteratiw Analysis, Prentice-Hall, Englewood Cliffs, N. J. (1%2).
15. H. Greenspan, et a]., Computing Methodr in Reactor Physics, Gordon and Breach. New
1
($f +f(r)s(r) = d3rs t(r)+(r).
v
(5-3 19) York (1968).
16. G. E. Forsythe and C. B. Moler, Computer Solution of Linear Algebraic Systems,
Prentice-Hall, Englewood Cliffs, N. J. (1967).
Since this must hold for a n y choice of S ( r ) a n d St(r), we will use it to our 17. T. Craig, SLODOG, A Modified One-group, One-dimensional Diffusion Code, Univer-
advantage by specifying S ( r ) as a unit point source a t ro: sity of Michigan Nuclear Engineering Report (1968).
18. E. Wachspress, Iteratioe Solution of Elliptic System, Prentice-Hall, Englewood Cliffs, N.
J. (1966), p. 22.
19. Ibid., p. 24, 25.
a n d St(r) a s the cross section Zd(r) characterizing a n imagined detector placed in 20. R. G. Steinke, A Review of Direct and Iterative Strategies for Solving Multidimensional
the core. T h e n we find Finite Difference Problems, University of Michigan Nuclear Engineering Report (1971).
21. W. R. Cadwell, et al., PDQ-An IBM-704 Code to Solve the Two-dimensional Few-
group Neutron Diffusion Equations, WAPD-TM-70 (1957) and related reports, WAPD-
TM-179, WAPD-TM-230, WAPD-TM-364, WAPD-TM-678.
22. D. L. Delp, et al., FLARE, A Three Dimensional Boiling Water Reactor Simulator,
Hence in this instance the adjoint flux is simply the response of a detector in the General Electric Company Report, GEAP 4598 (1964).
core t o a unit point source inserted a t a position r,. Once again we find that +t(r,J 23. E. G. Adensam, et al., Computer Methods for Utility Reactor Physics Analysis, Reactor
is a measure of the importance of a neutron event (in this case, the production, and Fuel Processing Technology, Vol. 12, No. 2 (1%9).
rather than the absorption, of a neutron) at a point ro in contributing t o the 24. E. Wachspress, Iterative Solution of Elliptic Systems, Prentice-Hall, Englewood Cliffs, N.
response of a detector with cross section Xd(r) (as opposed to reactivity). J. (1966). p. 83.
25. R. Froehlich, in Mathematical Models and Computational Techniques for Analysis of
F o r the simple one-speed diffusion model we have been studying, the adjoint
Nuclear Systems, USAEC CONF-730414-P2 (1973), p. VII-I.
flux +t(r) is identical t o the flux itself. Hence perturbations affecting the creation o r
destruction of neutrons will have the most pronounced affect i n those regions in
which the flux is largest. This is not true, however, for the more general multigroup PROBLEMS
diffusion model, a s we will see in Chapter 7.
5-1 Compare the derivation of the one-speed neutron diffusion equation with that for the
equation of thermal conduction, taking care to point out the assumptions and
REFERENCES approximations used in each case. Refer to any text on heat transfer such as those
listed at the end of Chapter 12.
1. J. W. Dettman, Mathemalical Methoak in Physics and Engineering, McGraw-Hill, New 5-2 By considering a plane source or absorber of neutrons located at the origin of an
York (1969). infinite medium, derive the interface condition Eq. (5-15) on the neutron current
2. B. Friedman, Principles and Techniques of Applied Mathematics, Wiley, New York (1956). density by modeling the source term in the one-dimensional diffusion equation as
3. 0.Arfken, Mathemorical Methods for Physicists, 2nd Edition, Academic, New York S&(x)and then integrating this equation over an infinitesimal region about origin.
(1970). 5-3 Compute the rms distance ( ( x ~ ) ) " ~a neutron will travel from a plane source to
4. P. M. Morse and H. Feshbach, Methods of Theorerical Physics, Vol. I, McGraw-Hill, absorption using one-speed diffusion theory. Compare this result with the rms
New York (1953), p. 173.
distance to absorption in a strongly absorbing medium (in which neutron scattering
5. Ibid., pp; 115-1 17.
can be neglected). In particular, plot the rms distance to absorption in water in which
6. B. Davison, Neutron Tramport Theory, Oxford U. P. (1958), pp. 51-55.
boron has been dissolved against the boron concentration to determine whether the
228 / THE ONESPEED DIFFUSION MODEL OF A NUCLEAR REACTOR THE ONESPEED DIFFUSION THEORY MODEL / 7.29
diffusion theory result for ((x2))'/' ever approaches the result characterizing a purely where a denotes the surface of the region. Yet another useful parameter characterizing
absorbing medium. (Use the thermal cross section data in Appendix A.) interfaces is the ratio of the current density J to the flux at the interface
It is possible to derive an expression for the relaxation or diffusion length L from the
one-speed transport equation characterizing a homogeneous medium
Determine the relation between these parameters and the albedo, assuming that
(we have assumed isotropic scattering for convenience). Seek a solution of the form diffusion theory can be used to describe the material adjacent to the region of interest.
v(x,p)=~(p)exp(-x/L) to this equation in order to eliminate the x-dependence. (It should be remarked that one frequently uses these concepts to characterize very
The resulting homogeneous equation for ~ ( p can ) be reduced to an algebraic
highly absorbing regions such as fuel elements or control rods in which diffusion
equation for L by eliminating J ?ldpX(p). By following this procedure, derive a
theory will usually not be valid in the highly absorbing region.)
transcendental equation for the diffusion length L: In reactor analysis it is frequently of interest to determine the neutron flux in a
so-called unit juel ceN of the reactor, that is, a fuel element surrounded by a
moderator. As a model of such a cell characterizing a cylindrical fuel element,
consider a fuel pin of radius a surrounded by a moderator of thickness b . For reasons
that will become more apparent in Chapter 10, one assumes that the fission neutrons
that slow down to thermal energies appear as a source uniformly distributed over the
Using the assumption that Z,CX,, expand L as a power series in X,/Lt, substitute moderator-but not directly in the fuel. Furthermore it is assumed that there is no net
this expansion into the equation above, and evaluate the coefficients of the expansion transfer of neutrons from cell to cell-that is, the neutron current vanishes on the
in order to derive Eq. (5-32) and obtain the transport corrections to the diffusion boundary of the cell (although the neutron flux will not vanish there).
length L -(D/L,)'/'. Determine the neutron flux in this cell geometry. In particular, determine the
Determine the neutron flux in a sphere of nonmultiplying material of radius R if an thermal utilization characterizing the cell by computing the fraction of those neutrons
isotropic point source of strength Soneutrons per second is placed at the center of the slowing down that is absorbed in the fuel.
sphere. Assume the sphere is surrounded by a vacuum. Consider a one-dimensional slab model of a fuel cell in which the center region
The Milne problem: Imagine a diffusing medium in the half space x > 0 with a source consists of the fuel, and the outer regions consist of a moderating material in which
of infinite magnitude at infinity such that the boundary condition on the flux is that neutrons slow down to yield an effective uniformly distributed source of thermal
+(x)-Soexp(x/L) as x+m. Perform the following calculations: neutrons Soneutrons/cm3~sec. Determine the neutron flux in this cell. In particular,
compute the so-called selfshielding factor j, defined as the ratio between the average
(a) Using one-speed diffusion theory and the boundary condition of zero reentrant
flux in the fuel to the average flux in the cell.
current, determine the flux in the medium.
Consider two isotropic point sources located a distance a apart in an infinite
(b) Repeat the solution of this problem using the extrapolated boundary concept.
nonmultiplying medium. Determine the neutron flux and current density at any point
(c) Determine the conditions under which these two boundary conditions might be
expected to yield similar results. in a plane midway between the two sources.
Determine the infinite medium Green's functions or diffusion kernels characterizing
Consider a slab of nonmultiplying material with a plane source at its origin emitting
So neutrons/cm2.sec. By solving this problem first with the condition of zero- cylindrical and spherical geometries.
reentrant current and then extrapolated boundaries, compare the absorption rate in By representing a plane source as a superposition of isotropic point sources, construct
the slab predicted by these two approaches. Also calculate the rate at which neutrons the plane source kernel GP,(x,x1) by using the point source kernel Gpt(r,r').
leak from the slab in each case. Obtain an expression for the plane source diffusion kernel characterizing a finite slab
Consider a thermal neutron incident on a slab-shaped shield of concrete I m in of width a by solving for the neutron flux resulting from a unit plane source at a
thickness, and determine the probability that: (a) the neutron will pass through the position x' in the slab. This can be most easily accomplished by seeking a separate
solution on either side of the source plane which satisfies the vacuum boundary
shield without a collision, (b) it will ultimately diffuse through the shield, and (c) it
conditions at either end, and then matching these solutions at the source plane using
will be reflected back from the shield. (For convenience, treat the concrete as if it had
the interface condition of continuity of flux and a discontinuity in the current density
the composition of 10% H20, 50% calcium, and 40% silicon.)
given by the source strength.
Consider an infinite nonmultiplying medium containing a uniformly distributed
Consider a neutron source emitting a monodirectional beam of neutrons into an
neutron source. If one inserts an infinitesimally thin sheet of absorber at the origin,
infinite medium. Using one-speed diffusion theory, calculate the neutron flux in the
determine the neutron flux throughout the medium.
medium. For convenience, locate your coordinate system with its origin at the source
Derive the expression given by Eq. (5-47) for the albedo characterizing a slab of
and align the x-axis along the source beam. Since the source is highly anisotropic, you
material of thickness a. In particular plot this albedo for a slab of water for various
cannot apply diffusion theory directly. Rather, compute the distribution of first
thicknesses. (Use thermal cross section data.) Comment on the behavior of the albedo
scattering collisions of the source neutrons along the x-axis, and then assume that
as given by Eq. (5-47) for both very thin and very thick slabs.
each of these collisions acts in effect as an isotropic point source of neutrons for the
One defines the blackness coeflicient characterizing a region as
subsequent diffusion theory analysis (assuming that such scattering is isotropic).
Use the method of variation of constants to determine the flux in a finite slab that
contains a uniformly distributed neutron source.
Construct the spatial eigenfunctions of the Helmholtz equation in spherical geometry.
230 / THE ONESPEED DIFFUSION MODEL & A NUCLEAR REACTOR
5-22 Construct the spatial eigenfunctions of the Helmholtz equation in infinite cylindrical
gwmetry.
1 THE ONESPEED DIFFUSION THEORY MODEL
Determine the thickness a for criticality when: (a) inner region is vacuum and (b)
/ 231
inner region is a medium with k, = I and same D and L as the outer slabs.
5-23 Construct the spatial eigenfunctions of the Helmholtz equation characterizing a A bare spherical reactor is made with 235U uniformly dispersed in graphite @= 1.7)
parallelepiped geometry. with an atomic ratio N,-/N2s= 1v. For the cross section values given below, calculate
5-24 Demonstrate explicitly that the eigenfunctions for the slab geometry are indeed the critical size and mass of the reactor according to one-group diffusion theory. If the
orthogonal.
reactor is modified by placing a cavity (vacuum) of half the total radius in its center,
5-25 Consider a slab of nonmultiplying material containing a uniformly distributed find the critical size for this case. Recalculate the critical radius if the center void is
neutron source. Determine the neutron flux in the slab: (a) by directly solving the
filled instead by a perfect absorber. Use as data: aF=4.3 b, aF=0.003 b, 0:~=105 b,
diffusion equation and @) using eigenfunction expansions. Demonstrate that these
a? 5 584 b, and D = 0.9 cm.
solutions are indeed equivalent.
A bare spherical reactor is to be constructed of a homogeneous mixture of D 2 0 and
5-26 Determine the neutron flux in a long parallelepiped column of nonmultiplying
235U.The composition is such that for every uranium atom there are 2000 heavy water
material caused by a uniform source distributed across the face of the column. In
molecules (i.e., ND10/N2s=2000). Calculate: (a) the critical radius of the reactor
particular determine the spatial behavior of the flux far from the source plane.
Describe how one might measure the neutron diffusion length in the column by using one-speed diffusion theory (Data: t12s=2.06, DDp-0.87 cm, XF0=3.3x
studying the spatial behavior of the flux. cm-', a?O=O.OO1 b, and ot5=678 b.) and (b) the mean number of scattering
collisions made by a neutron during its lifetime in this reactor.
5-27 The objective of this problem is to write a computer code to calculate the flux in a
uniform nonmultiplying slab containing an arbitrary source distribution. Perform the There is strong motivation to obtain as flat a power distribution as possible in a
following steps: reactor core. One manner in which this mag be accomplished is to load a reactor with
a nonuniform fuel enrichment. To model such a scheme, consider a bare, critical slab
(a) Derive the finite difference form of the appropriate one-speed diffusion equation
reactor as described by one-speed diffusion theory. Determine the fuel distribution
using vacuum boundary conditions. Use 30 mesh points and equal mesh spacing.
N d x ) which will yield a flat power distribution P(x)=wfZXx)cMx)=constant. For
(b) Write the equations in matrix form, _A9= g.
convenience, assume that fuel only absorbs neutrons and that it does not significantly
(c) Derive the steps necessary to solve this equation using the simultaneous relaxation scatter them. Also assume that all other materials in the core are uniformly distri-
method. (Gaussian elimination would be better, but less instructive.) buted.
(d) Write the necessary computer program, Input should consist of slab thickness, D, A one-dimensional slab reactor system consists of three regions: vacuum for x <O; a
Z,, and S,. Output should include a tabulation of xi, 4, and S,. Pay particular multiplying core for O< x < a; and an infinite nonmultiplying reflector for x > a.
attention to your convergence criterion and initial flux guess +i", for the inner Calculate the core thickness a that will yield a critical system.
iterations of the SR method. Suppose a cannot be made large enough to achieve criticality. Then determine the
(e) "Check out" your code by solving at least two problems with known analytical flux at all points when an external source So is uniformly distributed throughout the
solutions (e.g., cosine source or uniform source). Show how the inner iterations of reflector region x > a.
your program converge by plotting +,(") for various values of n. Consider a bare slab reactor with material composition such that 2,s 0.066 cm-',
5-28 Repeat the analysis of the timedependent slab reactor given in Section 5-3, but with 030.90, and vZ,=0.070 cm-'. Modify the one-speed diffusion computer code
the addition of a plane source of neutrons of strength Solocated at the origin of the developed in Problem 5-27 so that you can calculate the width a that will yield
slab. In particular determine the long time behavior of the reactor when it is critical. criticality.
5-29 Determine the geometric buckling B: and critical flux profile in the following bare One possible procedure is to guess an initial slab width a and then perform a source
reactor geometries: (a) sphere, @) infinite cylinder, and (c) parallelepiped. iteration calculation to determine ker To simplify the calculation, choose S(O)(x)= 1.0
5-30 We can define the power-peaking factor for a given reactor core as the ratio between =constant. The integrated fission source that appears in the estimate of k(") given by
the maximum power density and the average power density in the core. Recognizing Eq. (5-274) can be performed using simple trapezoidal quadrature
that the power density is proportional to the neutron flux in the one-speed approxima-
tion, compute the power-peaking factor for three common geometries: (a) sphere, (b)
cube, and (c) finite cylinder.
5-31 A homogeneous one-speed bare reactor has a cylindrical configuration. Determine:
(a) the radius and length of the reactor as functions of the buckling so that the volume
of the critical reactor, and hence its mass, is a minimum and (b) the minimum volume
as a function of buckling. where Ax is the mesh spacing. Again use 30 mesh points and require a convergence
5-32 Jezebel is a bare, fast, spherically shaped critical reactor constructed of pure 2 3 9 P ~ criterion on k(") of
metal (density 15.4 g/cm3). Calculate the critical radius and critical mass of the
reactor using the one-group data: v = 2.98, of= 1.85 b, a, ~ 0 . 2 6b, and a,,= 6.8 b.
5-33 There has been considerable interest in the possibility of super heavy nuclei with mass
numbers A >300. Such nuclei would be characterized by large values of 9 (-6-10).
Using the results of Problem 5-32, study the effect of varying v on the critical radius
After each criticality calculation, readjust the slab width a and recalculate key..After
and mass of a bare sphere of such material. several such calculations, plot kef, against a to determine the critical slab width a*.
5-34 Two infinite slabs, each of thickness a in the x direction, are separated by an inner Compare this with the analytical expression for a*.
region of thickness 2a and are bounded by vacuum on their outer surfaces. The slab
5-40 Investigate the convergence of the inner and outer iterations in the one-speed
material is of composition to give k,= 1.2 and thermal diffusion length of 50 cm.
diffusion code developed in Problem 5-39 for the following modifications:
232 / THE ONESPEED DIFFUSION MODEL OF A NUCLEAR REACTOR
(a) Use the outer iteration source S(")(x) as the initial estimate in the inner iterations.
(b) Determine the sensitivity of both inner and outer iterations to the convergence
criteria.
(c) Attempt to accelerate the outer iterations by source extrapolation. Nuclear Reactor Kinetics
5-41 Verify the general first-order perturbation theory result Eq. (5-312).
5-42 Why is the adjoint system introduced in developing the perturbation equations?
Illustrate your answer with an example showing that only the use of the adjoint
system will yield the desired result.
543 Describe a reasonable experimental procedure by which one could measure the
variation of neutron importance within a reactor core.
544 Calculate the error in the critical mass of a bare homogeneous spherical reactor due to
a 1% error in k. Assume that only the core size would be adjusted to give criticality.
5-45 Calculate the relative worth of a control rod bank inserted axially into a cylindrical
reactor core.
546 Consider a critical bare slab reactor of thickness a which is composed of a homo-
geneous mixture of fuel and moderator. Estimate the reactivity change if a thickness
8x of the fuel-moderator mixture at a position x is replaced by pure moderator. What
6x at a distance from the centerline of x-0.4 a is required to give the same reactivity
change as a perturbation thickness 6x0 at the center of the slab?
547 Variational methods can be used in a manner very similar to perturbation theory to
estimate the multiplication of a given core configuration using only crude guesses of
the flux in the core. For example, a useful variational expression for the multiplication
of a core described by a one-speed diffusion theory is
For a nuclear reactor to operate at a constant power level, the rate of neutron
Compare the accuracy of such a scheme for a slab of width a where the estimates of production via fission reactions should be exactly balanced by neutron loss via
the flux or "trial functions" +(x) are taken as simple quadratic polynomials that absorption and leakage. Any deviation from this balance condition will result in a
vanish on the boundaries of the slab [i.e., +(x)= I - ( 2 ~ / a ) ~ ] . time-dependence of the neutron population and hence the power level of the
reactor. This may occur for a number of reasons. For example, the reactor operator
might desire to change the reactor power level by temporarily altering core
multiplication via control rod adjustment. Or there may be longer term changes in
core multiplication due to fuel depletion and isotopic buildup. More dramatic
changes in multiplication might be caused by unforeseen accident situations, such
as the failure of a primary coolant pump or a blocked coolant flow channel or the
accidental ejection of a control rod.
It is important that one be able to predict the time behavior of the neutron
population in a reactor core induced by changes in reactor multiplication. Such a
topic is known as nuclear reactor kinetics. However, we should recognize that the
core multiplication is never completely under the control of the reactor operator.
Indeed since multiplication will depend on the core composition, it will also
depend on other variables not directly accessible to control such as the fuel
temperature or coolant density distribution throughout the reactor, but these
variables depend, in turn, on the reactor power level and hence the neutron flux
itself. The study of the time-dependence of the related processes involved in
determining the core multiplication as a function of the power level of the reactor
is known as nuclear reactor dynamics and usually involves a detailed modeling of
the entire nuclear steam supply system. Although we briefly discuss several of the
more important "feedback" mechanisms involved in determining core multiplica-
tion later in this chapter, our dominant concern is with predicting the time
behavior of the neutron flux in the reactor for a given change in multiplication.
233
234 / THE ONESPEED DIFFUSION MODEL OF A NUCLEAR REACTOR NUCLEAR REACTOR KINETICS / 235
The principal applications of such an analysis are not only to the study of I. THE POINT REACTOR KINETICS MODEL
operating transients in reactors, but also to the prediction of the consequences of
accidents involving changes in core multiplication, and to the interpretation of A. The Importance of Delayed Neutrons in Reactor Kinetics
experimental techniques measuring reactor parameters by inducing time-dependent
changes in the neutron flux. One can roughly distinguish between two different In our earlier treatment of the simple bare slab reactor (Section 5-111-C), we
types of analysis depending on the time scale characterizing changes in the neutron found that the neutron flux in such a system could be written as a superposition of
population or core properties. For example, one is interested in relatively short- spatial modes (or eigenfunctions) characteristic of the reactor geometry, each
term changes possibly ranging from fractions of a second up to minutes in length weighted with an exponentially varying time-dependence:
when analyzing normal changes in reactor power level (e.g., startup or shutdown)
or in an accident analysis. By way of contrast, changes in core composition due to
fuel burnup or isotope buildup usually occur over periods of days or months.
Needless to say, the analysis required for each class of time behavior is quite Here the spatial eigenfunctions were determined as the solution to the eigenvalue
different. problem [Eq.(5-232)l:
The reader will recall that we have already considered a particularly simple
example of nuclear reactor kinetics when we discussed the time behavior of the
neutron flux in a slab reactor model in Chapter 5. Although this earlier analysis
was useful for deriving the condition for reactor criticality, it is not valid for an
accurate description of nuclear reactor kinetics, since it assumed that all fission while the time eigenvalues A,, were given by
neutrons appeared promptly at the instant of fission. As we demonstrate in the next
section, it is essential that one explicitly account for the time delay associated with
delayed neutrons in describing nuclear reactor kinetics.
Fortunately the one-speed diffusion model we have been using to study reactor These eigenvalues are ordered as - A , > -A2> . . Hence for long times the flux
criticality is also capable of describing qualitatively the time behavior of a nuclear approaches an asymptotic form
reactor, provided we include the effects of delayed neutrons. Indeed such a model
is frequently too detailed for practical implementation in reactor kinetics analysis
due to excessive computatidn requirements, particularly when the effects- of
~henomenasuch as temverature feedback are included. For that reason. we will where we identify
begin our study of nuclear reactor kinetics by reducing the one-speed diffusion
model still further with the assumption that the spatial dependence of the flux in
the reactor can be described by a single spatial mode shape (the fundamental I =[oXa(l+ L ~ B ; ) ] - ' = ~ ~ ~lifetime
II of neutron in reactor,
mode). Under this assumption, we can remove the spatial dependence of the
vX,/Za k,
diffusion model and arrive at a description involving only ordinary differential k =- = ---=multiplication factor.
equations in time. This model is sometimes known as the point reactor kinetics I + L~B; I +L~B;
model, although this is somewhat of a misnomer since the model does not really
treat the reactor as a point but rather merely assumes that the spatial flux shape It would be natural to inquire as to just how long one would have to wait until such
does not change with time. Although we will rely heavily on the point reactor asymptotic behavior sets in. We can determine this rather easily by assuming that
kinetics model in our study of nuclear reactor time behavior, we will indicate its the reactor is operating in a critical state such that A,-0, and then estimating A,. If
generalizations to include a nontrivial spatial dependence as well as feedback we recall that for a slab B:= n2 ( n / q 2 , then
mechanisms.
There are two other aspects of nuclear reactor kinetics that we will not touch
upon in this chapter. The first topic involves the study of nuclear reactor control
and includes not only the analysis of the various schemes used to adjust core Now in a typical thermal reactor, a'-300 cm, c-3x l d cmjsec, and D-1 cm.
multiplication but those mechanisms by which changes in the core power level can Hence the higher order A,, are of the order of IW1000 sec-I, which implies that
affect multiplication as well. The second topic involves the study of long-term the higher order spatial modes die out very rapidly indeed.
changes in the core power distribution due to fuel depletion and the buildup of We can utilize this fact to bypass our earlier separation of variables solution of
highly absorbing fission products in the reactor core. the one-speed diffusion equation [Eq. (5-187)] by merely assuming a space-time
However the study of both of these subjects involves only a steady-state analysis
separable flux of the form
of the neutron flux in the reactor--or, at most, a sequence of steady-state criticality
calculations. Only the time-dependence of the slowly varying changes in core
composition such as those due to fuel depletion must be explicitly considered.
Therefore such topics can be more appropriately developed in later chapters. where +,(r) is the fundamental mode or eigenfunction of the Helmholtz equation,
2.36 / THE ONE-SPEED DIFFUSION MODEL OF A NUCLEAR REAmOR NUCLEAR REACTOR KINETICS / 237
Eq. (6-2). If we substitute this form into the one-speed diffusion equation, we find characterizing all fission neutrons as
that n(t) satisfies
In this sense, n(r) can be interpreted as the total number of neutrons in the reactor
at time t [if we normalize jd3r+,(r)= I]. Actually since the normalization of n(t) is Using the delayed neutron data given in Table 2-3, we find that this average
arbitrary, we could also scale the dependent variable n(r) to represent the total lifetime is typically (I)-0.1 sec, which is considerably longer than the prompt
instantaneous power P(r) being generated in the reactor core at any particular neutron lifetime 1-10-~- sec. Hence delayed neutrons substantially increase
time. That is, we could let the time constant of a reactor so that effective control is possible.
This fact suggests a related idea; suppose we consider a reactor that is very
slightly subcritical when only prompt neutrons are considered. Suppose further that
the fraction 8 of delayed neutrons provides just enough extra multiplication to
where w, is the usable energy released per fission event. Since the reactor power achieve criticality. This fraction will, in fact, control the criticality-and hence the
level is usually a more convenient variable to monitor, we will frequently express time constant. However if k - I > 8 , the reactor will be critical (or supercritical) on
the reactor time-dependence in terms of P(r). prompt neutrons alone, and the reactor period should become very short, since the
Equation (6-8) represents a somewhat simplified form of a more general set of delayed neutrons are not needed to sustain the chain reaction. Obviously we should
equations, which we shall derive in a moment, that are commonly referred to as the design a reactor such that this situation will never occur.
point reactor kinetics equations. This term arises since we have separated out the
spatial dependence by assuming a time-independent spatial flux shape +,(r). In this B. Derivation of the Point Reactor Kinetics Equations
sense we have derived a "lumped parameter" description of the reactor in which
the neutron flux time behavior is of the form of the product of a shape jacror +,(r) In actual fact, one cannot proceed so heuristically. We must first set up a set
and a time-dependent amplitude factor n(t), of equations describing the time dependence of the delayed neutrons. To this end,
we must define the precursor atomic number density:
Ci (r, t)d3r r expected number of "fictitious" precursors of
ith kind in d3r about r that always decay by (6- 13)
characterized by a time constant
emitting a delayed neutron.
T-E I period.
-reactor = (6-1 1 ) Note that Ci(r, t) is only some fraction of the true precursor isotope concentration,
since only a fraction of the ith isotope nuclei eventually decays by delayed neutron
This model is of course identical to that developed in our qualitative discussion of emission. For example, the ' l ~ rprecursor described in Section 2-11 is characterized
chai~i-reactionkinetics in Section 3-I-B except that one-speed diffusion theory has by a fictitious precursor concentration
now given us an explicit expression for the neutron lifetime I [Eq. (6-5)].
However there is, of course, something very important missing from this model. C,,(r,r) = (.029) (+7)Br(r, r ) . (6-14)
For by assuming a fission term in Eq. (5-186) of the form vZ&(r,f), we have
implied that all fission neutrons appear promptly at the time of fission. But we One can immediately write down a balance relation for these precursor concentra-
know that a very small fraction (8-0.7%) of such neutrons are emitted with tions by referring to our earlier discussion of radioactive decay to identify
appreciable time delay (recall the discussion of Section 2-11). Although these
I
delayed neutrons are only of minor significance in steady-state critical reactors,
I
Number of precursors
they are extremely important for reactor time behavior. For if we recall that the decaying in =A,Ci (r,t)d3r, (6- 15)
prompt neutron lifetime I is typically of the order sec in a thermal reactor
(10-' sec in a fast reactor), then it is apparent that the reactor period predicted by d3r/sec
I
this model would be far too small for effective reactor control.
I
Number of precursors
We can give a crude estimate of the influence of delayed neutrons on the reactor
time behavior by noting that the effective lifetime of such neutrons is given by their being produced = /3ivX&(r,t)d3r. (6- 16)
prompt neutron lifetime I plus the additional delay time A'; characterizing the in d 'r/sec
8-decay of their precursor. If we weight both prompt and delayed neutrons by their
respective yield fractions, we can estimate an average neutron lifetime (I) [This latter relation assumes the precursors don't migrate or diffuse before decay-
238 / THE ONESPEED DIFFUSION MODEL OF A NUCLEAR REACTOR NUCLEAR REACTOR KINETICS / 239
ing.] Hence our precursor concentration balance equation is just deviation of core multiplication from its critical value k = 1,
aci
-=
at -hiCi (r, t) + Pivz&(r, t). (6- 17)
Now our old friend, the one-speed diffusion equation, can still be used to describe
the flux-provided we treat the delayed neutron contribution to the fission source Notice here that we have explicitly indicated that k and hence p may be functions
explicitly by writing it as of time. We will comment more on this in a moment.
These definitions allow us to rewrite, Eqs. (6-21) in perhaps their most conven-
tional form
Hence our system of equations describing the neutron flux in a reactor including
delayed neutrons is
usually appear with energies below the fast fission threshold. We can take account
of these effects by characterizing each delayed neutron group by a slightly different important, one is usually forced to solve the time-dependent neutron diffusion
fast nonleakage probability and a resonance escape probability, say PLNLand p '. equation directly (at considerable expense). We will return to discuss such spatially
Then we would merely modify the definition of the parameters appearing in the dependent kinetics problems later in this chapter.
Of course one could proceed in a much more empirical fashion by noting that
point reactor kinetics equations as
the point reactor kinetics equations hold for more general situations than we have
considered, but then simply to postulate that the correct values for the parameters
8, A, and p are available (perhaps from experimental measurement). In this sense,
all detailed considerations of the flux spatial shape are avoided (or, rather, swept
under the carpet).
Let us now briefly examine the reactivity p(t)=[k(t)- I]/k(t) appearing in the
point reactor kinetics equations. We know that the multiplication factor k, and
hence p, depends on the size and composition of the reactor. In our specific
one-speed diffusion model of a bare reactor core, we have found an explicit form
for k [Eq. (6-5)]. Hence by changing the size or composition-say by inserting or
withdrawing a rod of absorbing material or adjusting a poison concentration-we
can change p and hence control the reactor. In this sense, p will in general be a
In lar& thermal power reactors, 6 is typically several percent larger than the function of time partly under the control of the reactor operator.
physical value 8, (although in small research reactors, may be as much as However for any reactor operating at power, p will also depend on the flux itself
20-30% larger than a). It should also be kept in mind that the delayed neutron due to several factors. First, the power level will influence the temperature of the
components of the reactor core. However the atomic concentrations of materials in
yield fractions /.Ii must be calculated as suitable weighted averages over the relevant
mixture of fuel isotopes. In fact one will find that in most thermal spectrum the core depend sensitively on their temperature. As the temperature changes, they
reactors these averaged delayed neutron fractions will decrease with core life since may contract or expand or change phase. This in turn will cause a change in the
the most common bred materials, ')'Pu and 2 3 3 ~ ,have somewhat lower delayed macroscopic cross sections-and hence in the reactivity. Furthermore temperature
neutron yields than ''U. changes may directly affect the microscopic cross sections (e.g., via the Doppler
In 239~u/2MU-fueledfast reactors, the situation becomes much more compli- effect). Finally, the atomic concentration of materials in the core will vary as
cated since there are as many as six fissioning isotopes, each characterized by six fission products are produced or fuel nuclei are fissioned and depleted. This will
different precursor groups. To handle 36 delayed neutron precursor concentration also strbngly influence the reactivity.
equations would be very complicated indeed. Fortunately the recent trend toward Such processes whereby the reactor operating conditions will affect the criticality
direct experimental measurement of isotopic yields as opposed to precursor group of the core are known as feedback effect and play an extremely important role in
data should allow one to eventually solve directly for the isotopic concentrations of reactor operation. Stated mathematically, such feedback effects imply that the
the major precursor isotopes (e.g., "Br, "Br, and ')'I) and only lump the minor reactivity must be regarded as a nonlinear function of the power level p[n(t),t].
isotopes into effective precursor groups. This would then eliminate the problem of Hence the point reactor kinetics equations are actually a coupled set of nonlinear
isotopic averaging for several fissionable isotopes, and decouple the data used in ordinary differential equations that are extremely difficult to analyze analytically,
the kinetics calcGation more effectively from the particular reactor core composi- with the exception of certain very simplified model cases.
tion. Although we will eventually discuss the physics of several of the more prominent
Perhaps the most serious approximation involved in the point reactor kinetics feedback mechanisms, we will initially limit our study to those cases in which p(t)
equations involves the assumption that the flux can be adequately represented by a is a specified function of time (and hence the point reactor kinetics equations are
single, time-independent spatial mode $,(r). It should first be noted that this shape linear)., Such a situation is commonly referred to as the zero-power point reactor
. function is actually not the fundamental mode characterizing a critical system, but model, since it ignores the feedback that would occur due to variations in the
rather the fundamental mode characterizing the reactor core that has been sub- reactor power level. Although limited in this sense, this model does reveal quite a
jected to a reactivity change away from critical. Nevertheless it is common to bit about the time behavior of the neutron population in a reactor.
utilize a shape function $,(r) characterizing a critical core configuration if the
reactor is close to a critical state or on a truly asymptotic period.
When the changes in core composition are sufficiently slow, as in fuel depletion 11. SOLUTION OF THE POINT REACTOR KINETICS
or fission-product poisoning studies, one can perform an instantaneous steady-state EQUATIONS
criticality calculation of the shape function $,(r), even though this shape will slowly
change with time. Such a scheme is known as the adiabatic appr~ximation.~ A. Solution With One Effective Delayed Group
More elaborate procedures exist for including a time variation in the shape We will first apply the point reactor kinetics equations to the situation in
function? However for rapidly varying transients in which spatial effects are which the reactivity is a prescribed function of time. In fact we will begin with the
242 / THE ONE-SPEED DIFFUSION MODEL OF A NUCLEAR REACTOR NUCLEAR REACTOR KINETICS / 243
very simple situation in which we imagine a reactor operating at some given power This set of homogeneous equations has a solution if and only if
level Po prior to a time, say t=O, at which point the reactivity is changed to a
nonzero value p,. Rather than attempting to solve the full set of point reactor
kinetics equations for this situation, we will first consider the case in which all
delayed neutrons are represented by one effective delayed group, characterized by
a yield fraction
The point reactor kinetics equations for this simplified case become and hence our general solutions will be of the form
To determine the unknown coefficients we can apply both the initial conditions
and the equations (6-35) to obtain four algebraic equations for the four unknowns
Note that we have chosen to work with the instantaneous reactor power level P(t) PI, P,, C,, and C,. Since these results are still rather cumbersome, we will examine
-
as our dependent variable. The precursor concentration is then slightly modified to
C = C,, = wfZfCgId.NOWprior to t 0, the reactor is operating at a steady-state
power level P,,. Hence we find that for t < O we must require
the solution in the case in which (/3-po+A~)2>4AXpo(or hA/P<l) and Ip0(<P.
Then the two roots are approximately given by
We have sketched this solution in Figure 6-1 for both positive and negative
This simple initial value problem can be solved in a variety of ways (the easiest reactivity insertion of an amount Jp01=0.0025 into a reactor characterized by
being by Laplace transforms as discussed in Appendix G). We will use a more P=0.0075, A-0.08 sec- I , and A = lo-' sec. In particular, it should be noted that
pedestrian approach by merely seeking exponential solutions of the form after a very rapid initial transient (s;'-.2 sec), the reactor time response becomes
exponential with a period of T = s;'-25 sec. This time constant characterizing the
P (t) = Pes', C (t) = Ce", (6-34) more slowly varying asymptotic behavior is occasionally referred to as the stable
reactor period.
where P, C, and s are to be determined. Then if we substitute these forms into Eq.
(6-31), we find the algebraic equations
B. The Inhour Equation
The reciprocal time constants s, and s, characterizing the single delayed
group model were given as the roots of a quadratic equation (6-36) in terms of p,,
p, A, and A. We can generalize this result to the situation in which there are several
groups of delayed neutrons by first rewriting Eq. (6-36) in a slightly different form.
244 / THE ONE-SPEED DIFFUSION MODEL OF A NUCLEAR REACTOR NUCLEAR REACTOR KINETICS / 245
Using our definitions of p and A in Eqs. (6-22) and (6-23), we can rewrite Eq. It should be noted that the range of p, is bounded by unity,
(6-36) as
This last limit is particularly interesting because it implies that no matter how much
negative reactivity we introduce, we cannot shut the reactor down any faster than
This equation, which expresses the decay constants s/ as the roots of a seventh- on a period T = 1 / 7 1 , determined by the longest-lived delayed neutron precursor. In
order polynomial, is known in reactor theory as the inhour equation. [This 235Ufueled thermal reactors, A;'-80 sec.
terminology arises from the very early attempts to express reactivity in units of Several other limiting cases are of particular interest:
"inverse hours" or inhours, which were defined as the amount of reactivity required
to make the reactor period equal to one hour. Reactivity is more commonly (1) Small reactivity insertions (po<P): Then we can assume that the magni-
expressed today either in decimals or percentages or pcm (per cent mille- lo-') of tude of so is small such that
Ak/k. We will later introduce another unit called the dollar, in which p is measured
in units of the delayed neutron fraction j3.1
The roots of Eq. (6-42) are most conveniently studied using graphical techniques.
In Figure 6-2 we have plotted the right-hand side of Eq. (6-42) for various values of Hence in the inhour equation (6-42) we can neglect sl in comparison with
s. The intersection of these curves with the line p=po yields the seven decay I - I and A, to write
constants 5 characterizing the time-behavior
equation (6-42) as
Asymptotic period (T), sec
$ 10-I
3
a
(noting that D(7)dr is just the probability that a delayed neutron will be emitted in (It will cross our path again.) Notice in particular that the reactivity insertion that
a time dr following a fission event at r=0) then we arrive at an "integro- gives rise to a purely sinusoidal power variation is periodic-but not sinusoidal (at
differential" form of the point reactor kinetics equations least for large power variations). One can show, in fact, that p(t) has a negative
bias. This fact proves of some importance in reactor oscillator experiments in
which an absorbing material is oscillated within the core and the core power
response is then measured in an effort to infer system parameters (e.g., P and A).
We can now rearrange this equation to yield p(t) in terms of P(t): 2. REACTIVITY AFTER A POSITIVE POWER TRANSIENT
As a second application of Eq. (6-52), suppose we determine the reactivity
present following a power transient in which the reactor power increases from its
initial value Po at I =0, and then decreases back to Po at a later time to
This particular relationship is important for two reasons: (a) it can be used in
principle to determine the time-dependence of the applied reactivity required to
yield a specific power variation-that is, to program the control rod motion and (b)
the interpretation of the measured power responses in transient analyses of reactiv- Now we know that for a positive transient, P(to-T)> P(to)= Po and hence the
ity changes can be used to provide information about the feedback mechanism in integral is always positive. Furthermore it is apparent that at t=to, the slope
the reactor. Several specific applications of this relationship are of considerable d P / d ~ l <, ~0. Thus we can infer that the reactivity following a positive power
interest: excursion must, in fact, be negative. That is, the reactor must be taken subcritical
to return the power to its original level.
I . PERIODIC POWER VARIATION
Suppose we seek the reactivity that will yield a sinusoidally varying power of 3. RAMP REACTIVITY INSERTION
the form Suppose the reactor power level is found to increase very rapidly from an
initial level Po as p0exp(at2) for t>0. Then using Eq. (6-52) we can find that for
long times I,
(Note that we must require Po> P, since a negative power would offend our
physical intuition.) Inserting this expression into Eq. (6-52) and performing a few
manipulations (left lor the reader's enjoyment to the exercises at the end of the Hence p(r) approaches a linear function of time-that is, a ramp insertion of
chapter), we arrive at reactivity above prompt critical. We can turn this calculation around to infer that
the response of the reactor power P(t) to a ramp insertion p(t)= yt should behave
as e ~ p [ ~ t ~ / 2for
A long
] times. This conclusion has particular relevance for certain
classes of fast reactor accident models in which fuel melting is assumed to lead to
core reassembly in a prompt supercritical configuration with a ramp reactivity
insertion. Obviously the power increase resulting from such an occurrence would
be very rapid indeed.
where D. Approximate Solutions
+=arg{ Y (iw)-'1. (6-55) As we have seen, exact solutions of the point reactor kinetics equations are
known for only a few special reactivity insertions. Hence we now turn our attention
and Y(iw) is a function which looks suspiciously like a chunk of the inhour to approximate schemes for solving these equations in the absence of feedback.
250 / THE ONESPEED DIFFUSION MODEL OF A NUCLEAR REACTOR NUCLEAR REAmOR KINETICS / 251
1. CONSTANT DELAYED NEUTRON PRODUCTION RATE one in which the prompt neutron lifetime is essentially taken to be zero such that
APPROXIMATION the power level jumps instantaneously to its asymptotic behavior (recall Figure
In certain problems, such as when the reactor is shut down by rapid insertion 6-1). This so-called prompt jump approximation is effected by merely neglecting the
of the safety rods, we are interested in the response of the reactor power to a given time derivative dP/dt in the point reactor kinetics equation
reactivity insertion in short time intervals following a time r,. During such short
time intervals we can ignore the change in the rate of production of delayed
neutrons, replacing C,(t) by Ci(0). Hence in this approximation, the point reactor
kinetics equation becomes
Since the delayed neutron production cannot respond immediately to a step change
where the effective source of delayed neutrons Q(t) is now presumed known and in reactivity, this model predicts that a reactivity jump from p, to p, causes an
given b, instantaneous change in reactor power from P I to P, as given by
After the scram rods have been fully inserted, the negative reactivity becomes
constant, and one can determine P(t) for subsequent time using our earlier For example, for a ramp insertion p(t)= ypf, the prompt jump approximation
solutions of the point reactor kinetics equation for constant p(t)= -p,. implies
2. THE PROMPT JUMP APPROXIMATION
In our earlier study of the response of the reactor power P(r) to a step change
in reactivity, we found that there was initially a very rapid transient behavior on a The prompt jump approximation is frequently found to yield an adequate
time-scale characteristic of the prompt neutron lifetime, followed by a more slowly description of reactor kinetic behavior. Numerical solutions of the point reactor
varying response governed by delayed neutron behavior. Since this transient is so kinetics equations have demonstrated the approximation is valid to within about
rapid, a very useful approximation to make for systems below prompt critical is I% up to reactivity insertions of p,=$0.50.6
252 / THE ONESPEED DIFFUSION MODEL OF A NUCLEAR REACXOR NUCLEAR REACTOR KINETICS / 253
then substitution of Eq. (6-70) into Eq. (6-69) and use of the identity
Hence to determine the behavior of the reactor power, we need only study the
poles of Z (s) and xs). However there is a great deal more we can do by employing
the very powerful methods of linear systems analysis. In particular we can study
yields (after a bit of manipulation) the stability of the reactor when it is operating at power (i.e., when we introduce
feedback).
The concept of a transfer function is much more general. The response (or
output) of any physical system to a signal (or input) applied to it can be expressed
in terms of such transfer functions. More precisely, we define
Thus far our analysis has been exact, merely consisting of a bit of algebraic
manipulation. We now introduce our key approximation by assuming that p(t) and Laplace transform of response -
p(r) are sufficiently small that we can neglect the second-order term p(t) p(r) to (6-77)
Transfer function = Laplace transform of input = Z (s).
write
case in which we allow feedback effects: As an example of the use of the reactor transfer function, let us determine the
response of the reactor to a sinusoidal reactivity input
~ ( r=) 6 ~ s i or.
n (6-8 1)
Output
Then if we note that the Laplace transform of this input is just
Feedback
' n p u i ~
where +(a) = arg[Z (io)] is known as the phase shijr of the transfer function, while
IZ(io)l is known as the system gain. The first term in Eq. (6-83) arises from the
poles of c(s) on the imaginary axis at s = ? io, while the remaining terms are due to
the poles of Z(s). These latter poles are just the roots sj of the inhour equation
Y (sj) = 0. Of course for the critical system we are considering, s, < s, < . . < s, =O.
Hence as t+m, only the oscillating terms and the s,=O term remain, and we find
the asymptotic behavior of the power oscillations as
Here, %(I) is the so-called unir impulse response, that is, the power response to a
unit impulse reactivity insertion. The 5 are the roots of the inhour equation, while
we have noted that for a critical system, s,,=O.
One can now use the convolution theorem (Appendix G) to reexpress Eq. (6-75)
in the "time domain" as
Notice in particular that there is a shift in the average power level of PoSp/oA.
If a reactor operating at a steady-state power level is subjected to a sinusoidal
perturbation in reactivity, the power will oscillate with the source frequency, but
with a phase shift +(o)=arg(Z) (actually a phase lag) and an amplitude propor-
Note that % ( I - 7 ) is just the Green's function for the linearized point reactor
tional to G (w) = ( Z(iw)l. Hence we can obtain the values of Z (s) on the imaginary
kinetics equation. Hence it is not surprising that Z(r) plays an extremely important
axis in the complex s-plane by measuring experimentally the amplitude and relative
role in the study of nuclear reactor kinetics. In particular one can infer the stability
phase of power oscillations as functions of frequency, induced by a sinusoidal
properties of a linear system by studying Z(r).
reactivity insertion with constant amplitude. This is the basis of the zero-power
One defines such a linear system to be stable if its response to any bounded input
pile-oscillator experiment. We have sketched the gain and phase shift of the zero
is also bounded. It is possible to show (see Problem 6-27) that a necessary and
power transfer function for a 235U/Z38U-fueled system in Figure 6-4.
sufficient condition for stability is for
models of the remainder of the plant suffices for on-line system performance
studies in the plant itself (using hybrid digital-analog process computers).
The most stringent application of the point reactor kinetics equations is in the
analysis of hypothetical reactor accident situations in which reactivity insertions
may be quite large (several S) and power transients very rapid. Indeed one
frequently encounters situations in which such equations are not valid. Typically
one uses them to provide a parameter study of the effects of various input data
such as reactivity insertions on reactor behavior under the postulated accident
conditions. Occasionally spatially dependent kinetics calculations are performed to
check the predictions of the point reactor model. Obviously an effort is always
made to generate conservative predictions.
For most severe transients in thermal reactors, such as a loss of coolant accident
or the ejection of a control rod, the usual transient thermal-hydraulics models are
adequate to supplement the core neutronics calculations. However in fast systems
very severe transients may lead to fuel melting and slumping with an appreciable
change in core configuration. Indeed such core rearrangements may produce large
positive reactivities, leading to very large energy release, followed by core disas-
sembly (that is, a small nuclear explosion). Such severe accidents, although quite
remote in probability, must nevertheless be studied in reactor analysis. They
require not only models of ,the core neutronics (e.g., the point reactor kinetics
equations) and thermal-hydraulics, but as well models of the physical motion of the
core due to melting and disassembly. Needless to say, such analyses can become
FIGURE W. Galn and p k sblft of tbe zero-power metor lmasfer f u a c t h very complicated indeed.
Most experimental measurements on reactors involving a reactivity change can 111. REACTIVITY FEEDBACK AND REACTOR DYNAMICS
be analyzed via the point reactor kinetics equations. Indeed since such experiments
are usually performed at zero power, it is frequently not necessary to append to the
point reactor kinetics equations suitable models of reactivity feedback. Although A. Mathematical Description of Feedback
we will briefly review a variety of experimental techniques commonly used to Thus far we have assumed that the reactivity p(r) appearing in the point
measure reactor kinetics parameters later in this chapter, we would merely mention reactor kinetics equations is a given function of time. However we know, in fact,
here as an example experiments in which the reactivity worth of control rods are that it depends on the neutron flux (or power level) itself. This dependence arises
measured by withdrawing a rod and then determining the resulting asymptotic because the reactivity depends on macroscopic cross sections, which themselves
reactor period. The inhour equation can then be used to relate this period to the involve the atomic number densities of materials in the core:
worth of the control rod.
Normal operating transients in nuclear power reactors are usually very easy to
analyze, since they are typically rather slow. This is because the rate of change of
reactor power level in a large plant is limited not by neutronic considerations, but
rather by the rate of temperature change allowable in nonnuclear plant com- Now it is easily understandable how the atomic density N(r, I ) can depend on the
ponents. For such slow transients (typically a 5% change in power per minute or reactor power level, since: (a) material densities depend on temperature T, which in
less for maneuvering) reactivity insertions are kept far below prompt critical, and
turn depends on the power distribution and hence the flux and (b) the concentra-
hence the prompt jump approximation, in which the mean generation time A is set tions of certain nuclei is constantly changing due to neutron interactions (buildup
equal to zero, is adequate. Furthermore on such a time scale the spatial power of poison or burnup of fuel). However it should also be noted that we have
shape usually remains in a fundamental mode and hence the point reactor kinetics explicitly written the microscopic cross sections as explicit functions of r and r .
equations are valid. The solution of the point reactor kinetics equations for such This must be done since the cross sections that appear in our one-speed diffusion
slow transients does necessitate consideration of reactivity feedback, however. model are actually averages of the true energy-dependent microscopic cross sec-
Frequently one finds the point reactor kinetics equations augmented with transient tions over an energy spectrum characterizing the neutrons in the reactor core. And
/ THE ONESPEED DIFFUSION MODEL OF A NUCLEAR REACTOR NUCLEAR REACTOR KINETICS / 259
this neutron energy spectrum will itself depend on the temperature distribution in
In this sense then, we define our incremental reactivities as
the core, and hence the reactor power level.
Such reactivity variation withtemperature is the principal feedback mechanism
determining the inherent stability of a nuclear reactor with respect to short-term
fluctuations in power level. Evidently our first task in constructing a model of
temperature feedback is to determine the temperature distribution in the core. Of
course, one could begin by writing the fundamental equations of heat and mass It is also useful to recall our definition of the incremental power
transport characterizing the reactor. That is, one could write the equations of
thermal conduction and convection using the fission energy deposition as the
source of heat generation. The motion of the coolant through the core would then From Eq. (6-88) we note that Spdp = 0] = SpJ Po]=0.
be described by the equations of hydrodynamics. As we will find in our detailed We can now sketch the block diagram characterizing a reactor with feedback as
development of thermal-hydraulic core analysis in Chapter 12, this fundamental shown in Figure 6-5. In particular, it should be recalled that we have already
approach results in a formidable set of nonlinear partial differential equations analyzed the "black box" describing the reactor kinetics, since this is just given by
unless further simplifications are introduced. the point reactor kinetics equation without feedback, which can be most con-
Hence it has become customary to avoid the complexity of a direct solution of veniently written in the form Eq. (6-72). This equation can be regarded as
the full set of thermal-hydraulic equations by replacing the spatially dependent
determining p(t) for any p(t). Of course we commonly analyze this equation under
description by a "lumped parameter" model similar in philosophy to the point
one of a variety of approximations, for example: (a) ignoring delayed neutrons
reactor kinetics model. The reactor core is typically characterized by several
(D rO)for very large reactivity insertions, (b) prompt jump approximation, and (c)
average temperatures, such as an average fuel temperature, moderator temperature,
linearization for small reactivity insertions, in which case we can solve Eq. (6-73) in
and coolant temperature. One then models a simple dependence of the reactivity
terms of the open loop or zero-power transfer function to find
on these temperature variables.
Of course one must append to the original point reactor kinetics model equations
describing the changes of these core temperatures produced by changes in the
reactor power level. Typically such calculations are based on a model of a single
average fuel element and its associated coolant channel, and describes the conduc- We therefore turn our attention to the study of the black box describing the
tion of fission heat out of the fuel element and into the coolant, and then its feedback functional Spdp].
subsequent convection out of the reactor core. We will consider such models in
some detail in Chapter 12. B. Models of Temperature Feedback
Our more immediate concern in this chapter, however, is to study how the results
of such thermal-hydraulic models, namely the core component average tempera- 1. TEMPERATURE COEFFICIENT OF REACTIVITY
tures such as T, (fuel) and TM (moderator or coolant), can be used in suitable
models of reactivity feedback. To this end let us return to consider our point Of most concern in the study of short-term reactivity feedback is the effect of
reactor kinetics model. We will write the reactivity p(t) as a sum of two contribu- the core temperature T on the multiplication of the core. One usually expresses this
tions in terms of a temperature coefficient of reactiuity a, defined as
The Sp notation signifies that the reactivity is measured with respect to the
equilibrium power level Po [for which p=O]. Furthermore, 6pex,(t) represents the
"externally" controlled reactivity insertion such as by adjusting a control rod. 7 -
r' '
This should more properly be referred to as an isothermal temperature coefficient, fission energy being released in the fuel, and hence the entire reactor core can
since it assumes that the core can be characterized by a single uniform temperature essentially be characterized by a single temperature. As the temperature of the core
T.
is increased, say by heating the primary coolant, there will be a decrease in core
Such a simple model illustrates quite clearly the effect of temperature feedback multiplication due to temperature feedback. One defines the so-called temperature
on reactor stability. For if a reactor were to possess a positive a,, then an increase defect of reactivity, Ap,, as the change in reactivity that occurs in taking the
in temperature would produce an increase in p, hence the power would increase, reactor core from the fuel-loading temperature (i.e., the ambient temperature) to
causing a further increase in temperature, and so on. In this sense, the reactor the zero power operating temperature
would be unstable with respect to temperature or power variations. The more
desirable situation is that in which a, is negatiw, since then an increase in T zero power aP
temperature will cause a decrease in p, hence a decrease in reactor power and AhD=/ -dT,
aT
temperature which tends to stabilize the reactor power level. T ambient
The temperature coefficient of reactivity depends on a great many processes
occurring within the core. Since the details of these processes depend on the The magnitude of APT, is primarily determined by the coolant temperature
specific reactor type under consideration, we will defer a discussion of just how coefficient of reactivity, and depends sensitively on the moderator-to-fuel ratio and
such temperature coefficients are evaluated to Chapter 14. soluble poison concentration. In a LWR ApTD-.02- .04(Ak/k).
However we should remark at this point that the concept of such an isothermal
temperature coefficient has a very limited range of usefulness in reactor analysis. In 2. POWER COEFFICIENTS OF REACTIVITY
a heterogeneous reactor. most of the fission energy release is confined to the fuel A far more useful parameter characterizing feedback is the power coefficient
elements. This energy must then be transferred by thermal conduction through the of reactivity defined by
fuel element, the clad, and then into the coolant before it can be withdrawn from
the core. Hence one has a very nonuniform temperature distribution. For example,
fuel centerline temperatures may range as high as 2700°C, while coolant tempera-
tures are as low as 250°C.
Hence one usually introduces several such temperature coefficients of reactivity,
where the T. are again the effective temperatures associated with each core
each of which characterizes the reactivity feedback due to a variation in the component. Such a parameter takes account of the temperature differences occur-
effective or average temperatures of each major component of the core, such as ring in a reactor while it is operating at power. The reactivity due to power
fuel, moderator, and structure feedback can then be written as
However one must also keep in mind that the various thermal processes in the core Obviously if a reactor is to be inherently stable against power-level fluctuations, it
are characterized by widely different time behaviors. The fuel temperature re- must be designed with ap<O.
sponds relatively rapidly to any power level changes. However it takes an appreci- A closely related quantity is the power defect, which is defined to be the change
able amount of time to transfer this energy to the coolant, and hence its tempera- in reactivity taking place between zero power and full power
ture response is much slower. For this reason, it is convenient to divide up the
temperature coefficient of reactivity into prompt and delayed components.
Effects that depend on the instantaneous state of the fuel-for instance, re-
sonance absorption (Doppler effect) or thermal distortion of fuel elements-may
be regarded as prompt, while effects that depend primarily on the moderator or The power defect can be quite sizable. For example, in the LWR the power defect
coolant-neutron energy spectrum and thermal expansion of moderator material- is typically of the order Ap,,+.Ol- .03(Ak/k).
are delayed. Prompt feedback mechanisms are of particular importance in reactor Thus far we have only discussed how the reactivity depends on temperature. The
safety studies, since they play a very important role in limiting any reactor remaining problem of how the temperature depends on the power level is strongly
transients that may occur. related to the reactor type and involves a thermal-hydraulics analysis of the core.
For these reasons, the isothermal temperature coefficient of reactivity is not For slow power changes, one can use the steady-state analysis developed in
really a very useful quantity. A more useful quantity is the change in reactivity Chapter 12 to determine the core temperatures for a given power level and hence
caused by a change in reactor power. We will consider this in detail in the next the power coefficient of reactivity in terms of the temperature coefficients of the
section. various components of the core.
There is one instance in which the concept of an isothermal temperature Such a steady-state thermal analysis will no longer be valid for more rapid power
coefficient is of some use, however. When the reactor is at zero power, there is no transients. For example, the thermal time constant of the fuel is frequently as large
262 / THE ONE-SPEED DIFFUSION MODEL OF A NUCLEAR R E A W R NUCLEAR REAmOW IUNFnCS / 263
as 10 seconds. For power transients on shorter time scales the fuel temperature response to an external reactivity insertion. First we will consider the power
coefficient of reactivity (principally determined by the Doppler effect) will be the oscillations resulting from a periodic reactivity insertion of small amplitude, since
dominant factor in determining the power coefficient of reactivity. this in effect measures the reactivity-to-power transfer function.
Although a realistic description of temperature feedback involves a transient A typical series of reactivity power transfer function gain measurement^'^^" is
thermal-hydraulics model of the reactor core such as those we will develop in shown in Fig. 6-6. In particular notice how different the transfer function at power
Chapter 12, several alternative lumped parameter models of the reactor tempera- is from the zero power transfer function. The marked resonance behavior in the
ture dependence are occasionally used for qualitative studies of reactor dynamics. vicinity of .03 cps at Po=550 kW in the first plot is apparent. Such behavior is due
One common model assumes a single effective coolant temperature T,, and to the presence of feedback. As the power level increases, the resonance peak
models the fuel temperature TFby Newton's law of cooling becomes narrower and higher. As we shall see, this implies that for sufficiently
large powers, the reactor is unstable.
Let us go back and develop a mathematical expression for the transfer function
characterizing a reactor with feedback.' We will restrict ourselves to small power
variations about the equilibrium level Po so that the feedback can be adequately
where K and Y are thermal constants characterizing the core. At the opposite
extreme would be an adiabatic model, in which the heat loss is assumed to be represented as a linear functional similar to Eq. (6-100)
negligible (such as in a very rapid transient)
Note that with this sign convention the negative feedback necessary for reactor
Still another model assumes consranr power remooal such that
All of these models can be considered as special cases of a general linear feedback
functional:
A little inspection should convince you that the feedback kernel h(t) in each of
these cases takes the form:
loo I I I I T
Newton's law of cooling: h(r) = aKe-7' 0.01 0.05 0.1 0.5 1
Adiabatic model: h(t)= aKIPo=O] (6-101) Frequency Icps)
stability will imply that h(r)<O. It should also be kept in mind that the feedback as in the zero power case, the inverse of L(s), /(I)= e-'{L(s)), is the Green's
kernel h(t) actually depends on the equilibrium power level P,. function for the point-reactor kinetics equation with feedback
If we substitute Eq. (6-102) into Eq. (6-72), we find
Let us examine L(s) in a bit more detail. Unlike Z(s), L(s) is analytic at s=O
with a value
Of course this equation is still nonlinear. We will consider only small power
variations p(r)< Po such that we can linearize Eq. (6-103) to write
Hence the long time response to a positive step reactivity insertion of magnitude
Spexl+Spo is just
p(rj= ~ ~ ~o ~ l b r 1 ( ~ ) 4 ~ ~ ~ p ~ (6-109)
~ ( 0 ) .
Now, as before, we will assume the reactor is operating at a steady state power
level Po prior to 1-0. Then by Laplace transforming Eq. (6-104) we find This implies that the reactor approaches a new equilibrium power level
This is in sharp contrast to the zero power (i.e., no feedback) reactor whose power
level grew exponentially for long times. The equilibrium state occurs when the
power level reaches a value such that the feedback reactivity just compensates for
the step reactivity insertion
where Z(s) is the usual zero power transfer function given by Eq. (6-76) while
H(s)- e ( h ( t ) ) is the jeedback transfer junction. We have further defined the
reacriuiry-to-power or closed-loop transfer function L(s) 2. RESPONSE TO A SINUSOIDAL REACTIVITY INSERTION
Let us once again consider a sinusoidal reactivity variation
&peXl(t)
= &posinwt. (6- 1 12)
Notice that as Po+O, L(s)-+Z (s), the zero power transfer function. This notation is
Then the long time response is given in analogy to the zero feedback case [Eq.
consistent with our earlier block diagram that has been relabeled in Figure 6-7. Just
(6-8411 by
'(')
-= 1( a )s i n (+ ) +(w)=arg(L(iw)]. (6-1 13)
Po
Now if we recall the form for L(s) given by Eq. (6-106), we can see that a
resonance in the gain G (w) = 1 L(iw)( will occur when
I - P,H (iw)Z (iw) = I + Pol H (io)Z (io)JexpiB+O, (6-114)
where
FIGURE 6-7.
u
Closed-loop Traoslw lunctbn 0 (a) = arg( - H (iw)Z (iw)).
7.66 / THE ONESPEED DIFFUSION MODEL OF A NUCLEAR REACTOR NUCLEAR REAaOR K I ~ E l W S / 267
However for the resonance condition [Eq. (6-122)) to be satisfied, we require that
negative feedback and hence shift the pole s, to the left and into the LHP.
both of the conditions
However the other poles will also shift. In particular, the pole s, will shift to the
right. For some sufficiently large power level P, the poles s, and s, will coalesce,
B (w) = 180' and POIH (iw)Z (iw)l= l
forming a complex conjugate pair that moves off of the real axis and into the
be simultaneously satisfied. These two conditions will determine a critical power complex plane. For still larger Po, this pair moves back to the right, until for some
level PC,, and a resonance frequency w, at which a true resonance situation will critical power level Po= PC,, the poles move into the RHP and the reactor becomes
arise. As long as the reactor power Po is kept below Pdl, only a finite resonance unstable. When the poles cross the imaginary axis, we observe a "resonance" in the
peak such as those illustrated in Figure 6-6 will occur. At this peak, the gain I L(iw)l transfer function gain IL(iw)J as we have already noted. Hence the problem of the
will assume a value linear stability of our equilibrium state reduces to the problem of determining the
sign of the real parts of the poles of L(s). In particular, we have found that a
reactor is linearly "strictly stable" when the poles of L(s) all have negative real
parts. The response of the power in the critical case, when any of these poles have
real parts equal to zero, is not correctly described by linear analysis and depends
In the light of this discussion we can more readily understand the transfer function on nonlinearities of the point reactor kinetics equation.
gain measurements given in Figure 6-6.
In almost all cases of reactor design, the critical power level PC,, above which the
reactor is susceptible to such instabilities is always quite far above actual or design
operating levels. Nevertheless it is important to be able to anticipate such inherent
I r-plane
- -l s-plane
in the complex s-plane. First notice that since Z(s) appears both in the numerator There are numerous tricks for determining the signs of these poles. These include
and denominator, its poles ,s [that is, the roots of the inhour equation Y(s,) methodsI2 such as those utilizing Nyquist diagrams, root-locus plots, and Routh-
= Z -'(s,)=O] "cancel." Hence the poles of L(s) are simply the zeros of Hurwitz criteria. However such subjects are more properly the concern of linear
system analysis or control theory, so we will simply refer the interested reader to
1- PoH (s)Z (s) = 0, the references listed at the end of the ~ h a p t e r . ' . ~ . ~
Now suppose that Eq. (6-119) were to have a simple root at s=so[i.e., a pole of
L(s)]. Then when we invert the Laplace transform, this pole will contribute a term 4. NONLINEAR POINT REACTOR KINETICS
in p(r) of the form exp(sot). Hence if so is in the right half s-plane (RHP), then p(t) Thus far our study of the point reactor kinetics equations with feedback has
will grow exponentially in time. This would imply an unstable response to an
been restricted to situations in which the reactivity changes and corresponding
applied reactivity perturbation (within the linear approximation, of course). If the power changes are sufficiently small that these equations can be linearized. In
root so lies in the LHP, then terms of the form exp(sot) will decay in time. Hence to particular our study of the stability of the reactor has been restricted to the
study reactor stability, it is obviously important to determine if any of the poles of
consideration of stability in the small-that is, to the study of perturbations and
L(s) [i.e., zeros of Eq. (6-1 19)] lie in the RHP.
responses sufficiently small for a linear analysis. However for larger perturbations
Instabilities may arise for sufficiently large power levels Po even with negative nonlinear effects must be taken into account. In these cases, the resulting con-
reactivity feedback, as the sequence of diagrams in Figure 6-8 indicate. For Po=O, clusions about stability may be quite different.
we can identify the poles of L(s) as just those of Z(s) [i.e., the inhour equation For example, we have seen that the linearized point reactor kinetics equation
roots]. Then, for the case in which a,= H(0)< 0, increasing Po will provide more predicts that the reactor will be unstable if the power exceeds some critical value
a / THE ONE-SPEED DIFFUSION MODEL OF A NUCLWR REACTOR NUCLEAR R E A O R KINETICS 1 269
P,,,. However even though the reactor is linearly unstable, it may be stable in the flux resulting from a source in a subcritical assembly is measured as fuel is added
nonlinear description. Hence it is of some interest to study the significance of linear to the assembly. If the reactor is characterized by a multiplication factor k, then
stability theory within the more general framework of nonlinear stability theory. one can crudely think of the amplification of the original source neutrons by the
It should be remarked, however, that there are many different approaches to assembly as given by
nonlinear point reactor kinetics-none of which is completely satisfactory.
Furthermore such results as do exist provide only sujficienr (as opposed to neces-
sary) conditions for stability. Sufficiency conditions are usually much too
restrictive for practical application. Furthermore, the very limited validity of the
models investigated by nonlinear stability methods usually discount their value for
practical reactor design. And of course there is also the more pragmatic philosophy where So is the rate at which source neutrons are emitted in the reactor. The usual
adopted in most nuclear reactor design favoring an ultraconservative design that procedure for a safe approach to delayed critical in core loading consists of
guarantees linear stability under all conceivable operating conditions, hence obviat- plotting M-' (or the reciprocal neutron counting rate) as a function of some
ing the need for a nonlinear stability analysis. parameter that controls reactivity (e.g., fuel mass) and then extrapolating this M - I
plot to zero to determine the critical loading (as sketched in Figure 6-9). During a
5. SOME FINAL COMMENTS ON REACTOR STABILITY ANALYSIS stepwise approach to delayed critical by the reciprocal multiplication method, the
It might seem that an integral part of any reactor safety analysis would be an neutron level following each addition of reactivity must be allowed to stabilize in
investigation of the stability of the reactor design. In particular such an analysis order to obtain an accurate indication of the asymptotic multiplication before
should consider the possibility that the reactor might become unstable under some proceeding with the next reactivity addition. As one approaches delayed critical,
feasible combination of operating conditions. In practice, however, such stability the buildup of the precursor concentrations can become bothersome.
studies do not play near as significant a role in reactor design as one might expect.
Reactor instabilities usually can occur only if one of the temperature feedback
coefficients happens to be positive over some range of operating conditions.
However it is usually quite easy to design a power reactor so that it will always be
characterized by large, negative temperature coefficients (as we will see in Chapter
14). Indeed all power reactor designs tend to be ultraconservative in this respect (as
they do with respect to all safety considerations). Hence there has been relatively
little motivation to perform detailed stability analyses of reactor core designs. And
in those instances in which methods of stability analysis have been applied to
practical reactor designs, investigations have usually been limited to the linear
domain.
Fuel mass FIGURE 6-9: Cdtial loading measurements
IV. E X P E R I M E N T A L D E T E R M I N A T I O N O F R E A C T O R
KINETIC PARAMETERS AND REACTIVITY
From our earlier development of the point reactor kinetics equations, we have 2. FUEL SUBSTITUTION TECHNIQUES'
found that the three most important parameters characterizing the kinetic behavior
of a nuclear reactor are the reactivity p, the prompt generation time A, and the One can determine reactivity by uniformly substituting a p .&sonfor a small
effective delayed neutron fraction 8. A variety of experimental techniques have fraction of the fuel. The poison is usually chosen such that it has the same
been developed to measure both these as well as dynamic (i.e., feedback) scattering and absorption properties as the fuel (or as closely as possible). Then
characteristics of the reactor. Such methods can be classified as either static or first-order perturbation theory yields a reactivity change due to the substituted
dynamic measurement techniques. poison of
B. Dynamic Techniques for Reactivity Measurements Using Eqs. (6-124) and (6-125), along with the equilibrium forms of the point-
reactor kinetics equation prior to the source jerk, one can find
1. ASYMPTOTIC PERIOD MEASUREMENTS
Perhaps the simplest type of kinetic measurement one can perform is to make
a small perturbation in the core composition pf a critical reactor, and then to
measure the stable or asymptotic period of the resultant reactor transient. Using
the inhour equation, one can then infer the reactivity "worth" of the perturbation
from a measurement of the asymptotic period. It should be noted that the period Hence the dollars subcritical of the assembly is given very simply by
method for all practical purposes applies only to positive periods, since negative
periods are dominated by the longest delayed neutron precursor decay and hence
provide very low sensitivity to negative reactivity.
2. ROD DROP METHOD The source jerk method is very similar to the rod drop method in concept. However
Let us consider a reactor operating at some equilibrium power level Po when it is somewhat simpler to perform since it requires only the removal of a small mass
it is suddenly shut down by the introduction of a negative reactivity - S p (the "rod of material (the source) in contrast to the rapid release of one or more control rods.
drop"). From our earlier studies of the point reactor kinetics equations we know
that after a few prompt neutron lifetimes, the reactor power level drops to a lower 4. ROD OSCILLATOR METHOD
level P , , determined by the amount of reactivity insertion and remains at this If one oscillates a control rod in a sinusoidal fashion in a critical reactor,
"quasistatic" level until it is ultimately decreased by delayed neutron precursor there will be a corresponding oscillation in the reactor power. By measuring the
decay. We can use our earlier result [Eq. (6-64)] from the prompt jump approxima- gain and magnitude of the phase shift characterizing the power oscillations, one
tion to write can measure the reactor transfer function, either at zero power or operating power:
Hence we can solve for the reactivity insertion in dollars in terms of the power For low powers, L(iw)+Z(iw). In this case, the high-frequency behavior of the
levels Po and PI as transfer function
To determine PI, one need only extrapolate back the asymptotic behavior follow- can be used to measure either the prompt generation time A or the reactivity worth
ing the rod drop to I =O. of the oscillating rod,
3. SOURCE JERK METHOD
A very similar technique can be used to measure the multiplication of a
subcritical assembly. Suppose we consider such a subcritical system maintained at
a power level Po by a neutron source of strength S,. One can express the Such oscillator measurements can also be used to determine the stability
equilibrium power level Po using the point reactor kinetics equation as characteristics of the reactor.
subject to the initial condition applying to the pulsed source, +(r,O)=+,(r). For In this case, the a, versus B: curve appears as that shown in Figure 6-10.
long times, the flux in the assembly will approach the asymptotic form One can use this technique to directly measure reactivity. For a rapid pulse, the
time behavior of the flux or power will be determined by prompt neutron kinetics
such as
where +,(r) is the fundamental mode o f . the assembly geometry. Hence the
asymptotic behavior of the flux is governed by the decay constant
Hence a, can be used to infer p, provided P and A are known. If the core is
maintained at delayed critical, then p = 0 and the pulsed neutron method measures
If one measures a, for various assembly sizes, then plotting a, against B: will yield
vZ, (the intercept at B;=0) and OD (the slope). (See Figure 6-10.) Actually, there
are higher order terms in B: due to both transport and energy-dependent correc-
tions
C. Noise Analysis in Nuclear R e a ~ t o r s ' ~ - ' ~
Thus far we have been analyzing the reactor as essentially a deterministic
which add a curvature to the a,(B;) plot. system. Yet we know that all of the dynamic variables describing the reactor
In multiplying assemblies if one assumes that both the pulse injection and (power level, flux, temperature, etc.) actually fluctuate in a statistical fashion about
measurement are performed on a time scale short compared to delayed neutron some mean value. We are unable to predict with certainty the future values of these
lifetimes, then one can use variables but rather can only specify the probability that they will assume a certain
value.
The statistical nature of the neutron diffusion process has been stressed re-
peatedly. Moreover the concept of a cross section, and, indeed, even the quantum
mechanical description of the neutron interactions, must be interpreted in a
to find decay constants of the form probabilistic sense. There are numerous other sources of such statistical fluctua-
tions or noise in a reactor, such as in fluid flow, coolant boiling, and mechanical
ao5 V(Z. - v 2 3 + VDB,~. (6- 136) vibration. In fact, statistics enters even into the measurement of the reactor state
due to detector noise.
At low power levels the statistical fluctuations associated with the fundamental
nuclear processes occurring in the core will be dominant. However at higher power
levels, the reactor noise d l be predominantly due to disturbances of a nonnuclear
nature (e.g., coolant flow). Regardless of its origin, reactor noise is important in
reactor dynamics for at least two reasons. First, it interferes with the precision with
which a desired quantity can be measured in a reactor; that is, one must extract the
signal of interest out of the background noise. Second, however, since the noise
originates from various processes occurring in the reactor, it can actually serve as a
source of information about the system. The latter of these properties is of most
interest to us here. We will study how noise analysis can be used to measure the
transfer function of the reactor.
If the functions x(t) and y(t) are periodic, then the limit process can be omitted then construct both the time-correlation functions ~ ( 1 and
) their Fourier trans-
provided T is chosen as the period. In general, however, x(t) and y(t) will not be forms numerically. For example, if measurements are taken at delay intervals of
periodic. In fact we will later consider them to be random wriables in the sense that AT, then the cross-spectral density is given by
only their probability distribution can be specified.
Now we will set up the cross correlation between reactivity and power (using a
notation in which the limit process T-tm is to be understood)
Hence by merely measuring the cross correlation, one can determine both the
amplitude and phase of the transfer function.
This experiment serves, then, as an alternative to a reactor oscillator transfer-
function measurement. Unlike the latter, it does not suffer from background noise
(rather taking advantage of such random fluctuations), and hence does not require
nearly so large an input signal. However both of these experiments suffer from the
fact that one must perturb the reactor by introducing an externally controlled
reactivity signal in order to perform the measurement. It is possible to bypass this
difficulty and measure the amplitude of L(io) directly from the inherent noise
where we have identified the autocorrelation function of reactivity naturally present in the reactor.
2. AUTOCORRELATION MEASUREMENTS
As a final step, we take the Fourier transform of Eq. (6-142), defined as Consider the autocorrelation of the fluctuations in the reactor power
to find
If we substitute in Eq. (6-107) for p(t)/P, and again take a Fourier transform in
time, we find
but S ( l ( u ) ) is just the closed loop transfer function L(iw). Hence we find that we
can write
Here ePp
= 9 ( q p p is) referred to as the cross-spectral density, while FPp= F(p$,) is Hence the square of the magnitude (i.e., the gain) of the transfer function can be
referred to as the reactivity (or input) spectral density. determined as the ratio of the power and reactivity spectral densities.
To apply this result in a practical measurement, one varies the reactivity of the Although we can easily measure the power spectral density, the reactivity
reactor in a random or pseudorandom fashion and then measures the correspond- spectral density characterizing inherent reactivity noise is not experimentally ac-
ing variations in the reactor power or flux. In this sense, the measurement is very cessible. If one has reason to believe that such fluctuations are truly random in
similar to that involved in a rod oscillator experiment. To construct the time- nature, then 9(9,) =constant, and the measurement of the power autocorrelation
correlation function, one measures both 6p,,,(r) and p(r) at the time r and at function by itself will yield the amplitude of the transfer function (although phase
.
various delayed times 1' + A T , t AT, .. . Using these measurements, one can
information will have been lost).
276 / THE ONE-SPEED DIFFUSION MODEL OF A NUCLEAR REACTOR NUCLEAR REACXOR KINETICS / 277
Such a measurement of IL(iw)l is extremely simple. It has the advantage that it time scale of the order of the effective neutron lifetime (I), the point reactor
does not perturb the reactor. Indeed it can be used in a real time mode to obtain an kinetics equation should be regarded as suspect.
instantaneous measurement of the reactor transfer function. Its only significant In these instances, one must take explicit account of the spatial dependence of
disadvantage is its limited accuracy. the neutron flux. In fact one is frequently forced to perform a brute force
Such noise measurements have been used to measure a large variety of nuclear numerical solution of the time-dependent neutron diffusion equation and precursor
reactor parameters, and have proven to constitute a valuable diagnostic technique equations. Unfortunately in most cases in which such spatial effects are significant,
in determining reactor kinetic behavior. one cannot rely on the one-speed approximation to provide an adequate descrip-
tion of the neutron energy-dependence. Hence a direct numerical study of nuclear
reactor kinetics1$usually involves the solution of the multigroup, time-dependent
V. SPATIAL EFFECTS IN REACTOR KINETICS diffusion equations-at a considerable computational expense.
For such finite difference solutions of the multigroup diffusion equations to be
feasible, it is necessary to employ a variety of numerical techniques to accelerate
Thus far our analysis of time-dependent nuclear reactor behavior has been
the computation. These calculations are usually coupled with transient thermal-