炉物理の研究 第71号 （2019年2月）

My two-cents experience in doing reactor physics
Yung-An Chao
November 2018 in Shanghai China

Professor Akio Yamamoto asked me if I could write an article for the Reactor Physics Division of
AESJ on my experience of doing reactor physics, mainly for the interest of students and young colleagues.
As I started to write it, I realized that I committed to a very challenging task when I accepted his invitation.
I decided to choose a few examples to convey my way of learning and doing reactor physics. As thetargeted readers are mainly students, I must apologize to more senior colleagues if the writing sounds
trivial or too much preaching.
1. Two examples in text book reading:First a very simple example. Resonance shielding is an important subject in reactor physics. To
explain why a sharply peaked absorption resonance will be less absorptive than a flatter and broader
resonance, text books typically show mathematically how a dipped flux spectrum results in lessabsorption via the evaluation of the reaction rate integral. I find it simpler to understand the physics
by thinking of using a roll of toilet paper to wipe out spilled water. It is obvious that pulling out the
paper from the roll will be much more effective than applying the roll directly. The outside of the roll“shields” the interior of the roll. The same analogy applies to the spatial shielding effect of a lumped
fuel rod embedded in moderator versus a dispersed mixture of fuel with the moderator.
Then a little more intriguing example. The core system eigenvalue (k) and the corresponding
reactivity are fundamental concepts in reactor physics. Text books tell us that a steadily operatingcore is at criticality, k=1. But a core cannot start by itself without some initial “extraneous” neutron
source (S). A simple point model for neutron balance gives that the neutron flux level is proportional
to S/(1-k), with k being the ratio of neutron production to absorption. Then how can the flux for a
critical core be maintained if the flux would become infinity at criticality? The explanation is that when
the core is brought close to criticality, the flux level increases rapidly to the targeted level while thecore still staying just a tiny bit subcritical, never really reaching criticality. Text books always present
an eigenvalue problem without including the S term in the equations, because when close to criticality
S is completely negligible compared to the chain reaction neutrons and its presence has nil effect on
the flux distribution. This simple point model also points to a few interesting things.
• A core steady state can occur only when it is sub-critical. When the core goes critical orsuper critical, it must be in a transient state. Otherwise the flux would become negative.
• If one writes a code to solve the diffusion equation in presence of S, one would run into
trouble if the core is very close to criticality. To avoid this problem, one could estimate
the eigenvalue on flight and switch to the eigenvalue problem when the S term becomes
negligible.
• Although the flux distribution near criticality is given by the fundamental mode
independent of the shape of S, in a deeply subcritical core the flux distribution may belargely affected by the S distribution. The deviation of the flux distribution from thefundamental mode distribution indicates the degree of subcriticality.
• Consequently, the effective neutron multiplicative capability of a subcritical core may bequite different from what the core eigenvalue would suggest. However, one should not
be misled that one might “gain” any reactivity safety margin because of this. The safety炉物理の研究 第71号 （2019年2月）

margin is dictated by how far the core is from criticality, which is independent of S anddepends solely on the eigenvalue.
2. On nodal methods:
Nodal method has been a dominating tool in core analysis over the past thirty plus years, andlikely to remain so for at least some years to come. Although I worked on it and used it myself, I havealways had some reservation about it. The most popular nodal method is the transverse integration
nodal method because of its simplicity and efficiency. The two pillars of the method are thediscontinuity factor and the transverse leakage profile, the former of which is semi-empirical while
the latter is quite ad hoc. Particularly the construction of the transverse leakage profile is quite
arbitrary, as theoretically there is even no reason why the profile should be continuous across theintersect of the four adjacent nodes. The assumption for the justification is that the transverse leakage
is supposed to be small, thus the result will not be so sensitive to its detailed profile. This may beadequate for diffusion calculation. But when using SP3, as to be discussed in the following section 5,
the correct interface and boundary conditions will contain higher order tangential derivatives, whichis the transverse leakage on the surface. In such a case if one would still use the ad hoc transverse
leakage profile, one would run into iteration instability and nonconvergence. A more rigorous solution
method must be used instead. As for the discontinuity factor, it typically comes from reflective unit
assembly homogenization, which is the weakest model and the biggest problem in today’s core
analysis method, and I believe it will eventually be the killer of today’s method.
Another tool typically used together with nodal method is CMFD (coarse mesh finite difference).
The classical finite difference current between two adjacent nodes is proportional to the flux
difference 𝜙𝐿 − 𝜙𝑅 . CMFD makes a simple revision to it by adding another term proportional to thesum 𝜙𝐿 + 𝜙𝑅 . The coefficient of the difference term is theoretically well known, while the coefficientof the sum term is empirically determined by brute force matching the CMFD neutron current to thecurrent from the higher level nodal method calculation. Although with no theoretical justification, thisrecipe works because of having one additional free parameter to play with. I did not feel comfortable
about it and hesitated to use it until I could understand it. I noticed that this recipe simply means ageneral linear combination of the two fluxes, i.e. 𝐷𝐿 𝜙𝐿 − 𝐷𝑅 𝜙𝑅 , and this form can be rigorously
derived in the 1D case with both DL and DR analytically determined. The argument goes as follows. In
1D case, the analytic flux solution in each of the two adjacent nodes has two free parameters. The
four free parameters in total for the two nodes are reduced to two by imposing the flux and current
continuity condition on the interface. Hence there is only one free parameter per node, which can bechosen as the node average flux. Therefore, it must be possible to analytically express the interface
current in terms of the two fluxes only. Pursuing along this line I derived the Analytic CMFD (ACMFD)
equations for the diffusion case in both Cartesian and hexagonal geometry. (I also derived the ACMFD
coefficients for the transport case, although I never had a chance to implement it.)
A popular use of CMFD in nodal method is the so called two nodes model, where the higher level
nodal calculation is done only for two adjacent nodes to determine the CMFD coefficients. The fullcore calculation is then done only with CMFD, which in turn provides the boundary value of the two
nodes problem. Personally, I think relying only on the (empirical) CMFD for the whole core calculation
is risky and I have a concern that the iteration may not be robust for all engineering problems. I prefer
to use CMFD/ACMFD purely as an acceleration tool such that thru the iteration process both the lower
level CMFD and the higher level nodal method calculation are done for the whole core problem. In

炉物理の研究 第71号 （2019年2月）

case the acceleration piece would cause any problem, the design code could automatically turn it off
to assure convergence.
3. On conformal mapping:
By late 1980’s the transverse integration nodal method had become the mainstream solver
because of its simplicity and efficiency in engineering applications. When people tried to apply themethod to hexagonal geometry, they ran into the difficulty of how to handle the resulting severe
singularity at the vertices of the hexagon. This problem caught my interest. I recalled what I read in
my school days that in classical electrostatics people used in early days conformal mapping to do thedesign of electric capacitor of different shapes because the Laplacian operator ∇2 in the electrostatic
equations is invariant under conformal mapping. Since diffusion equation has the same operator, weshould be able to use conformal mapping to map a hexagon to a rectangle to then apply the transverse
integration method. This should alleviate the singularity problem.
Since I did not expect to have the opportunity to work on this “academic” problem in ourcompany, I wrote a Letter to the Editor of Nuclear Science and Engineering to publicize this idea,
hoping that someone would pick it up. But nobody did, and I realized that I had to do more to savethis idea. I went to the company library to find a 1952 book “Dictionary of Conformal Representation”
and borrowed from a colleague a floppy disc of Fortran II to work out at my home the numericalcalculation of the mapping function. To be practical for implementation, the next step was toapproximate the numerical results with analytic functions in a closed form. At this point a manager
told my boss not to allow me wasting my time on academic work. Then an opportunity came just intime for me to make the work into a university cooperation project. To figure out the approximation
function, I used what I learned in school about complex variable function that the only function being
analytic everywhere (including at infinity) is a constant function. All other functions must havesingularities. Thus, a very powerful and efficient way of approximating a function is to locate and best
approximate its singularities.
Then the East European countries opened, and the company got interested in the refuelingbusiness of VVER reactors. Using the conformal mapping method, we could much more easily adapt
the existing core design codes to hexagonal cores. Interestingly the manager who did not like me
wasting time on academic work was then appointed a manager in the VVER project, and I believe hethen much appreciated the value of what I persisted on.
4. On Stiffness Confinement Method (SCM) in space-time kinetics:
Stiffness in kinetics is caused by the prompt neutron lifetime being orders of magnitude shorter
than the delayed neutron lifetime such that the time discretization has to take very small time steps
controlled by the prompt neutron lifetime. The idea of SCM was inspired by the observation that the
stiffness effect results in a prompt jump rise of the prompt neutron population while does not lead toany irregular behavior in the precursor response. This is because the precursor population is muchlarger while its population change is much smaller compared to the prompt neutron. By introducing
a dynamic frequency, the stiffness can be effectively confined only to the prompt neutron equation
which can be analytically solved after the stiffness absent precursor equation is solved. The idea was
first demonstrated and confirmed numerically with point kinetics.

炉物理の研究 第71号 （2019年2月）

When generalized to space time kinetics, the dynamic frequency is split into a shape term plus anamplitude term. The shape frequency is a function of both space and time, while the amplitude
frequency is a function of time only. The former is obtained from the prompt neutron shape change
while the latter from the prompt neutron total population change. The prompt neutron kinetic
equation is converted to an equivalent static eigenvalue equation, where the dynamic frequencies
divided by neutron velocity are treated as artificial macroscopic cross-sections. If the eigenvalue is notequal to unity, then an adjustment on the amplitude frequency is made, according to the resulting
point kinetic equation. This process is done iteratively until convergence. The flux distribution is thenupdated in both shape and amplitude using the corresponding dynamic frequencies. This method
allows the use of much larger time discretization steps while still giving accurate results. But unlike
other methods, when using SCM one needs to be cautious about using small time steps, which may
cause instability. The convergence criteria in the SCM eigenvalue problem calculation must becompatible with the actual change of the eigenvalue and eigen function between the two ends of thetime step. If the convergence residual is not negligible compared to the actual change in the time step,then the time step should be enlarged.
The above work was related to thesis research of my students. At the time our company did not
find any market incentive to convert the 3D kinetics code into a commercial product, because safety
analysis based on 3D kinetics was considered unnecessary. The code was put on shelf several years
without transient thermo-hydraulics added to it. It was unplanned and totally unexpected that its firstengineering application turned out to be the simulation of DRWM (see next section), where nothermo-hydraulics feedback is needed because the low power physics test is essentially done at zero
power. Sometimes one does not know what a research project would eventually lead to.
5. On Dynamic Rod Worth Measurement (DRWM):
The problem of DRWM development is how to use a point-model based reactivity meter toanalyze data from a 3D transient process. The key is of course how to account for the spatial effects.
A quick physics analysis concludes that there are two spatial effects, the static spatial effect and the
dynamic spatial effect. Consider a “Gedanken” experiment on a symmetric core with one control rod
completely inserted and its symmetric partner completely withdrawn. If we swap the positions of thetwo identical rods, there should be no reactivity change of the core, yet the excore reactivity meter
will record flux change and indicate positive or negative reactivity change pending on where the meter
is located. This is a static spatial effect, due to core configuration change per se. The dynamic spatial
effect arises because the delayed neutron distribution change trails behind that of prompt neutron.
Simulation calculation can predetermine the proper signal correction factors due to these two effects.
But then there is the subtle question whether the use of these correction factors would “mask” themeasurement data such that the measurement result would always agree better with the core designprediction. A qualitative physics analysis shows that this correlation is generally negative, which tends
to enlarge the discrepancy between the prediction and measurement. This anti-masking feature will
give more conservative results, which is acceptable from the safety point of view. A large amount ofsimulation calculation was done to confirm the conclusion of this qualitative analysis.
When we got the first DRWM testing data from a core startup physics measurement, thediscrepancy between the design prediction and measurement of rod worth was excessively large. The
management decided to abandon the development project. But I noticed that the problem came fromthe heavy rods, one set of which was unusually heavy with 1800pcm rod worth. A simple point kinetics
model indicates that the rod worth would be roughly proportional to the logarithm of the ratio of the炉物理の研究 第71号 （2019年2月）

initial flux to the final flux. When such heavy rods were inserted the signal dropped so low that thebackground could severely contaminate the final flux value and thus affect the results significantly.
Since the background signal should be the same for all the sets of rods, I best estimated a constantbackground to subtract from all the signals. All the discrepancy largely disappeared! The estimatedbackground was a fraction of the lowest signal, which was quite reasonable. The project survived, andthe background compensation has since been a crucial part of DRWM.
6. On the Generalized SPN theory (GSPN):
Despite its simplicity and practical usefulness, I was always a little skeptical about SPN because of
its lack of a theoretical foundation. Although I never used SP3 myself I started to think about it more
seriously when my student needed to calculate the SP3 discontinuity factor for cell homogenization.
To properly define the discontinuity factor, one needs to compare the reference transport solution to
the SP3 solution. But there is no way of deducing the reference value of the SP3 functions from thetransport solution because we don’t know how to relate the SP3 functions to the transport angular
flux, thus popping up the question “what is the angular flux representation corresponding to the SP3
solution?”
As I started to investigate this problem, I soon realized that the rigorous derivation of the SPN
differential equations for a homogeneous medium has never been a problem, which has been donein various ways by people. (To my big surprise, I discovered that the first derivation was as early as1944 by Mark and a little later by Davison when they developed the PN theory, although not calling it
by SPN at the time.) The problem has always been not knowing how to derive the interface andboundary conditions (IBC). In practice people have been using the (ad hoc) 1D type IBC on the surface
value and the normal direction derivative of the SPN solution functions. I then questioned why we caneasily derive the IBC for PN but not for SPN? Physically the IBC should follow from the continuity
requirement on the angular flux. In case of PN the angular flux is represented via spherical harmonics
expansion with the PN functions as the expansion coefficients. In such a case, angular flux continuity
is the same as the continuity of the PN functions. The derivation of IBC for PN therefore does not need
the reconstruction of the PN angular flux. For SPN, however, there is no such a simple expansion.
Therefore, in search for the SPN IBC we must go beyond the domain of the SPN functions to ask whatis the SPN angular flux representation corresponding to the SPN solution functions? I faced the samequestion again: we don’t understand the physical meaning of SPN unless we know its corresponding
angular flux representation.
I struggled with the quest for the SPN angular flux representation for a couple of years and it
eventually leaded me to the development of the GSPN theory, where I had to learn and use themathematical tools of solid harmonics, the Helmholtz decomposition of vector, the Davison lemma
for solid harmonics, and the generalized Helmholtz decomposition of tensor. (Ironically all these tools
were used or developed by Davison sixty years ago in his fascinating work of partially successfulattempt to solve the PN equations analytically.) The GSPN theory transforms and recasts the PN
equations into N+1 layers, where each layer contains a set of diffusion type differential equations. The
equations in each layer is decoupled from the other layers. Each layer contributes a piece of theangular flux which is explicitly given in terms of the spatial gradients of the layer’s solution functions.
The complete PN angular flux is the sum of all the pieces contributed by the different layers. Each time
when a layer is added, the summed angular flux changes such that the derived IBC changes as well.
Therefore, although the diffusion type differential equations in different layers are decoupled, the IBC
for the layers are nevertheless coupled depending on how many collected layers one wants to use.

炉物理の研究 第71号 （2019年2月）

The lowest layer in GSPN is identical to the traditional SPN except for a different IBC, which is no longer
1D type and contains higher order tangential derivatives in addition to normal derivatives. Since
conventional numerical methods have difficulty of handling the higher order tangential derivatives in
IBC, a new numerical method is developed to solve the GSPN equations.
Concluding remarks:
Usually the first thing I do when I start to work on a problem is to think about the physics behind sothat I can properly define the problem from the physics point of view. The crucial step is then to build acorresponding physics model, where necessary assumptions (i.e. approximations) are introduced.
Henceforth it can be converted to a mathematical problem, where I may have to find or learn newmathematical tools to solve it. I am not good at reading a lot of equations, where often I quickly get lost.
If I could first understand what the physics problem is, then I could more patiently focus on reading theequations repeatedly to understand it. If I cannot understand something in my own way, I havereservation on believing and accepting it.
A professor I most respected in graduate school told us “Physics is a science of approximation, andmathematics is a science of exactness. Mathematics is not physics.” When we solve a physics problem,
we always need to build a model with assumptions. Building the model is doing physics, the science ofapproximation. But following what my professor said, I would like to add “Science of approximation does
not mean approximate science.” Good physics should avoid ad hoc approximations.
Everyone has his/her own way of working and his/her own taste of appreciation. Undoubtedly myexperience is very subjective and may not even worth two cents to you. But anyway, I have fulfilled mycommitment to Prof. Yamamoto.
References:
CMFD/ACMFD:• Y. A. Chao, "A Theoretical Analysis of the Coarse Mesh Finite Difference Representation in
Advanced Nodal Methods", The Proceedings of M&C'99-Madrid, Mathematics and Computation,
Reactor Physics and Environmental Analysis in Nuclear Applications, Vol. 1, pp.117 - 126, Madrid,
Spain, September 27-30, 1999.
• Y. A. Chao, "Coarse Mesh Finite Difference Methods and Applications", Proceedings of the
PHYSOR 2000 ANS International Topical Meeting on Advances in Reactor Physics andMathematics and Computation into the Next Millennium (on CD), Pittsburgh, USA, May 7-12,
2000
Conformal mapping:• Y. A. Chao, "Comments on the Treatment of Transverse Leakage in Advanced Nodal Codesfor Hexagonal Nodes," Letter to Editor, Nucl. Sci. Eng., 1991, 109, p.423.
• Y. A. Chao and N. Tsoulfanidis, “Conformal Mapping and Hexagonal Nodal Methods – I:
Mathematical Foundation”, Nucl. Sci. Eng., 1995, 121, pp.202-209.
• Y. A. Chao and Y. A. Shatilla, “Conformal Mapping and Hexagonal Nodal Methods – II:
Implementation in the ANC-H Code”, Nucl. Sci. Eng., 1995, 121, pp.210-225.
Stiffness confinement method:

炉物理の研究 第71号 （2019年2月）••

Y. A. Chao and A. Attard, “A Resolution of the Stiffness Problem of Reactor Kinetics”, Nucl. Sci.
Eng., 1985, 90, pp.40-46.
Y. A. Chao and P. Huang, “Theory and Performance of the Fast-Running Multidimensional
Pressurized Water Reactor Kinetics Code, SPNOVA-K”, Nucl. Sci. Eng., 1989, 103, pp.415-419.

DRWM:• Y. A. Chao, D. M. Chapman, D. J. Hill and L. R. Grobmyer, “Dynamic Rod Worth Measurement”,
Nucl. Tech., 2000, 132, pp.403-412.
GSPN theory:• Y. A. Chao, “A new and rigorous SPN theory for piecewise homogeneous regions”, Annals of
Nuclear Energy, 2016, 96, pp.112-125.
• Y. A. Chao, “A new and rigorous SPN theory – Part II: Generalization to GSPN”, Annals of Nuclear
Energy, 2017, 110, pp.1176-1196.
• Y. A. Chao, “A new and rigorous SPN theory – Part III: A succinct summary of the GSPN theory, the(3)P3 equivalent GSP3 and implementation issues”, Annals of Nuclear Energy, 2018, 119, pp.310321.
• Y. A. Chao, L. Peng and C. Tang, “A new and rigorous SPN theory - Part IV: Numerical testing of(0)
GSP3 using the plane wave expansion method”, in preparation for submittal to Annals of
Nuclear Energy.