PyRK : A Python Package For Nuclear Reactor Kinetics

In this work, a new python package, PyRK (Python for Reactor Kinetics), is introduced. PyRK has been designed to simulate, in zero dimensions, the transient, coupled, thermal-hydraulics and neutronics of time-dependent behavior in nuclear reactors. PyRK is intended for analysis of many commonly studied transient scenarios including normal reactor startup and shutdown as well as abnormal scenarios including Beyond Design Basis Events (BDBEs) such as Accident Transients Without Scram (ATWS). For robustness, this package employs various tools within the scientific python ecosystem. For additional ease of use, it employs a reactor-agnostic, object-oriented data model, allowing nuclear engineers to rapidly prototype nuclear reactor control and safety systems in the context of their novel nuclear reactor designs.


Introduction
Time-dependent fluctuations in neutron population, fluid flow, and heat transfer are essential to understanding the performance and safety of a reactor.Such transients include normal reactor startup and shutdown as well as abnormal scenarios including Beyond Design Basis Events (BDBEs) such as Accident Transients Without Scram (ATWS).However, no open source tool currently exists for reactor transient analysis.To fill this gap, PyRK (Python for Reactor Kinetics) [Huff2015], a new python package for nuclear reactor kinetics, was created.PyRK is the first open source tool capable of: • time-dependent, • lumped parameter thermal-hydraulics, • coupled with neutron kinetics, • in 0-dimensions, • for nuclear reactor analysis, • of any reactor design, • in an object-oriented context.
As background, this paper will introduce necessary concepts for understanding the PyRK model and will describe the differential equations representing the coupled physics at hand.Next, the implementation of the data model, simulation framework, and numerical solution will be described.This discussion will include the use, in PyRK [Huff2015], of many parts of the scientific python software ecosystem such as NumPy [vanderWalt2011] for array manipulation, SciPy [Milman2011] for ODE and PDE solvers, nose [Pellerin2015] for testing, Pint [Grecco2014] for unit-checking, Sphinx [Brandl2009] for documentation, and Matplotlib [Hunter2007] for plotting.

Background
Fundamentally, nuclear reactor transient analyses must characterize the relationship between neutron population and temperature.These two characteristics are coupled together by reactivity, ρ , which characterizes the departure of the nuclear reactor from criticality: The reactor power is stable (critical) when the effective multiplication factor, k, equals 1.For this reason, in all power reactors, the scalar flux of neutrons determines the power.The reactor power, in turn, affects the temperature.Reactivity feedback then results due to the temperature dependence of geometry, material densities, the neutron spectrum, and reaction probabilities [Bell1970].This concept is captured in the feedback diagram in Figure 1.One common method for approaching these transient simulations is a zero-dimensional approximation which results in differential equations called the Point Reactor Kinetics Equations (PRKE).PyRK provides a simulation interface that drives the solution of these equations in a modular, reactor design agnostic manner.In particular, PyRK provides an object oriented data model for generically representing a nuclear reactor system and provides the capability to exchange solution methods from one simulation to another.
The Point Reactor Kinetics Equations can only be understood in the context of neutronics, thermal-hydraulics, reactivity, delayed neutrons, and reactor control.

Neutronics
The heat produced in a nuclear reactor is due to nuclear fission reactions.In a fission reaction, a neutron collides inelastically with Fig. 1: Reactivity feedback couples neutron kinetics and thermal hydraulics a 'fissionable' isotope, which subsequently splits.This reaction emits both heat and neutrons.When the emitted neutrons go on to collide with another isotope, this is called a nuclear chain reaction and is the basis of power production in a nuclear reactor.The study of the population, speed, direction, and energy spectrum of neutrons in a reactor as well as the related rate of fission at a particular moment is called neutronics or neutron transport.Neutronics simulations characterize the production and destruction of neutrons in a reactor and depend on many reactor material properties and component geometries (e.g., atomic densities and design configurations).

Thermal-Hydraulics
Reactor thermal hydraulics describes the mechanics of flow and heat in fluids present in the reactor core.As fluids are heated or cooled in a reactor core (e.g.due to changes in fission power) pressure, density, flow, and other parameters of the system respond accordingly.The fluid of interest in a nuclear reactor is typically the coolant.The hydraulic properties of this fluid depend primarily on its intrinsic properties and the characteristics of the cooling system.Thermal hydraulics is also concerned with the heat transfer between the various components of the reactor (e.g., heat generation in the reactor fuel heat removal by the coolant).Heat transfer behavior depends on everything from the moderator density and temperature to the neutron-driven power production in the fuel.

Reactivity
The two physics (neutronics and thermal-hydraulics) are coupled by the notion of reactivity, which is related to the probability of fission.The temperature and density of materials can increase or decrease this probability.Fission probability directly impacts the neutron production and destruction rates and therefore, the reactor power.The simplest form of the equations dictating this feedback are: The Point Reactor Kinetics Equations (PRKE) are the set of equations that capture neutronics and thermal hydraulics when the time-dependent variation of the neutron flux shape is neglected.That is, neutron population is captured as a scalar magnitude (a point) rather than a geometric distribution.In the PRKE, neutronics and thermal hydraulics are coupled primarily by reactivity, but have very different characteristic time scales, so the equations are quite stiff.
In the above matrix equation, the following variable definitions are used: β = fraction of neutrons that are delayed (8) β j = fraction of delayed neutrons from precursor group j (9) The PRKE in equation 5 can be solved in numerous ways, using either loose or tight coupling.Operator splitting, loosely coupled in time, is a stable technique that neglects higher order nonlinear terms in exchange for solution stability.Under this approach, the system can be split clearly into a neutronics sub-block and a thermal-hydraulics sub-block which can be solved independently at each time step, combined, and solved again for the next time step.

PyRK Implementation
Now that the premise of the problem is clear, the implementation of the package can be discussed.Fundamentally, PyRK is object oriented and modular.The important object classes in PyRK are: • SimInfo: Reads the input file, manages the solution matrix, Timer, and communication between neutronics and thermal hydraulics.• THComponent : Represents a single thermal volume, made of a single material, (usually a volume like "fuel" or "coolant" or "reflector" with thermal or reactivity feedback behavior distinct from other components in the system. • Material : A class for defining the intensive properties of a material (c p , ρ, k th ).Currently, subclasses include FLiBe, Graphite, Sodium, SFRMetal, and Kernel.
A reactor is made of objects, so an object-oriented data model provides the most intuitive user experience for describing a reactor system, its materials, thermal bodies, neutron populations, and their attributes.In PyRK, the system, comprised by those objects is built up by the user in the input file in an intuitive fashion.
Each of the classes that enable this object oriented model will be discussed in detail in this section.

SimInfo
PyRK has implemented a casual context manager pattern by encapsulating simulation information in a SimInfo object.This class keeps track of the neutronics system and its data, the thermal hydraulics system (THSystem) and its components (THComponents), as well as timing and other simulation-wide parameters.
In particular, the SimInfo object is responsible for capturing the information conveyed in the input file.The input file is a python file holding parameters specific to the reactor design and transient scenario.However, a more robust solution is anticipated for future versions of the code, relying on a json input file rather than python, for more robust validation options.
The current output is a plain text log of the input, runtime messages, and the solution matrix.The driver automatically generates a number of plots.However, a more robust solution is anticipated for v0.2, relying on an output database backend in hdf5, via the pytables package.

Neutronics
The neutronics object holds the first 1+j+k equations in the right hand side of the matrix equation in 5.In particular, it takes ownership of the vector of 1 + j + k independent variables and their solution.It also customizes the equations based on paramters noted in the user input file.The parameters customizing these equations for a particular reactor include α i for each component, j, Λ, k, and the fissionable nuclide.
The Neutronics class has three attributes that are sufficiently complex as to warrant their own classes: PrecursorData, Decay-Heat, and ReactivityInsertion.
A Neutronics object can own one PrecursorData object.In this class, the input parameters J and the fissionable nuclide are used to select, from a database supplied by PyRK, standardized data representing delayed neutron precursor concentrations and the effective decay constants of those precursors (λ d, j , β j , ζ j .That nuclear data is stored in the PrecursorData class, and is made available to the Neutronics class through a simple API.
A Neutronics object can also own one DecayHeat object.In this class, the input parameters K, and the fissionable nuclide are used to select, the fission product decay data (λ FP,k , ω k , κ k .The DecayHeat class provides a simple API for accessing those decay constants, fission product fractions, and weighting factors. Finally, a Neutronics object can own one ReactivityInsertion object.This defines the external reactivity, rho e xt, resulting from control rods, external neutron sources, etc.With this ReactivityInsertion object, the Neutronics class is equipped to drive a reactivity insertion accident scenario.That is, an accident scenario can be driven by an insertion of reactivity (e.g. the removal of a control rod).In PyRK, this reactivity insertion capability is captured in the ReactivityInsertion class, from which reactivity insertions can be selected and customized as in Figure 2.

THSystem
A reactor is made up of many material structures which, in addition to their neutronic response, vary in their temperature response.These structures may include fuel, cladding, coolant, reflectors, or other components.In PyRK, a heat transfer model of the changing temperatures and material properties of those components has been implemented as a lumped capacitance model.
This model approximates heat transfer into discrete components, approximating the effects of geometry for "lumps" of material.
In this model, heat transfer through a system of components is modeled analogously to current through a resistive circuit.Table 1 describes the various canonical forms of lumped capacitance heat transfer modes.
Based on the modes in Table 1, we can formulate a model for component temperatures specific to to the geometry of a particular reactor design.This might include fuel pellets, particles, or pebbles, cladding, coolant, reflectors or other structures in the design.
Fundamentally, to determine the temperature change in a thermal body of the reactor, we rely on relations between temperature, heat capacity, and thermal resistance.As in Table 1, the heat flow out of body i is the sum of surface heat flow by conduction, convection, radiation, and other mechanisms to each adjacent body, j [Lienhard2011]: Note also that the thermal energy storage and release in the body is accordingly related to the heat flow via capacitance: where Together, these form the equation:

THComponent
The THSystem class is made up of THComponent objects, linked together at runtime by heat transfer interfaces selected by the user in the input file: In the above example, the mat argument must include a Material object.

Material
The PyRK Material class allows for materials of any kind to be defined within the system.This class represents a generic material and daughter classes inheriting from the Material class describe specific types of material (water, graphite, uranium oxide, etc.).The attributes of a material object are intrinsic material properties (such as thermal conductivity, k t h) as well as material-specific behaviors.
Given these object classes, the burden of the user is then confined to:

Quality Assurance
For robustness, a number of tools were used to improve robustness and reproducibility in this package.These include: • GitHub : for version control hosting [GitHub2015] • Matplotlib : for plotting [Hunter2007] • Nose : for unit testing [Pellerin2015] • NumPy : for holding and manipulating arrays of floats [vanderWalt2011] • Pint : for dimensional analysis and unit conversions [Grecco2014] • SciPy : for ode solvers [Oliphant2007], [Milman2011] • Sphinx : for automated documentation [Brandl2009] • Travis-CI : for continuous integration [Travis2015] Together, these tools create a functional framework for distribution and reuse.

Unit Validation
Of particular note, the Pint package [Grecco2014]_ is used for keeping track of units, converting between them, and throwing errors when unit conversions are not sane.For example, in the code below, the user is able to initialize the material object with k th and c p in any valid unit for those quantities.Upon initialization of those member variables, the input values are converted to SI using Pint.

Minimal Example : SFR Reactivity Insertion
To demonstrate the use of this simulation framework, we give a minimal example.This example approximates a 1-second impulsereactivity insertion in a sodium cooled fast reactor.This type of simulation is common, as it represents the instantaneous removal and reinsertion of a control rod.The change in reactivity results in a slightly delayed change in power and corresponding increases in temperatures throughout the system.For simplicity, the heat exchanger outside of the reactor core is assumed to be perfectly efficient and the inlet coolant temperature is accordingly held constant throughout the transient.

Minimal Example Results
The results of this simulation are a set of plots, the creation and labelling of which are enabled by matplotlib.In the first of these plots, the transient, beginning at time t = 1s, is driven by a step reactivity insertion of 0.5 "dollars" of reactivity as in Figure 3.The power responsds accordingly as in Figure 4. Finally, the temperatures in the key components of the system follow the trends in Figure 5.
These are typical of the kinds of results nuclear engineers seek from this kind of analysis and can be quickly re-parameterized in the process of prototyping nuclear reactor designs.This particular simulation is not sufficiently detailed to represent a benchmark, as the effect of the cladding on heat transfer is neglected, as is the Doppler model controlling fuel temperature feedback.However, it presents a sufficiently interesting case to demonstrate the use of the PyRK tool.

Conclusions and Future Work
The PyRK library provides a modular simulation environment for a common and essential calculation in nuclear engineering.PyRK is the first freely distributed tool for neutron kinetics.By supplying and API for ANSI standard precursor data, a modular material definition framework, and coupled lumped parameter thermal hydraulics with zero-dimensional neutron kinetics in an object-oriented modeling paradigm, PyRK provides designagnostic toolkit for accident analysis potentially useful to all nuclear reactor designers and analysts.
where ρ(t) = total reactivity ρ f (t) = reactivity from feedback ρ ext (t) = external reactivity insertion and where , based on dT i dt and the external reactivity insertion.•THSystem : Manages various THComponents and facilitates their communication during the lumped parameter heat transfer calculation.

Fig. 2 :
Fig.2: The reactivity insertion that can drive the PyRK simulator can be selected and customized from three models.

•
defining the simulation information (such as duration or preferred solver) • defining the neutronic parameters associated with each thermal component • defining the materials of each component • identifying the thermal components • and connecting those components together by their dominant heat transfer mode.

Fig. 3 :
Fig. 3: A prompt reactivity insertion, with a duration of 1 second and a magnitude of 0.05δ k/k drives the simulation.It represents the prompt partial removal and reinsertion of a control rod.

Fig. 4 :
Fig.4: The power in the reactor closely follows the reactivity insertion, but is magnified as expected.

Fig. 5 :
Fig. 5: While the inlet temperature remains constant as a boundary condition, the temperatures of fuel and coolant respond to the reactivity insertion event.