Matrix elements from the spatial discretization.
A summary is described about nuclear power reactors analyses and simulations in the last decades with emphasis in recent developments for full 3D reactor core simulations using highly advanced computing techniques. The development of the computer code AZKIND is presented as a practical exercise. AZKIND is based on multi-group time dependent neutron diffusion theory. A space discretization is applied using the nodal finite element method RTN-0; for time discretization the ?-method is used. A high-performance computing (HPC) methodology was implemented to solve the linear algebraic system. The numerical solution of large matrix-vector systems for full 3D reactor cores is achieved with acceleration tools from the open-source PARALUTION library. This acceleration consists of threading thousands of arithmetic operations into GPUs. The acceleration is demonstrated for different nuclear fuel arrays giving extremely large matrices. To consider the thermal-hydraulic (TH) feedback, several strategies are nowadays implemented and under development. In AZKIND, a simplified coupling between the neutron kinetics (NK) model and TH model is implemented for reactor core simulations, for which the TH variables are used to update nuclear data (cross sections). Test cases have been documented in the literature and demonstrate the HPC capabilities in the field of nuclear reactors analysis.
- high-performance computing
- nodal finite element method
- parallel computing
- graphics processing unit
- neutronic-TH coupling
The mathematical models representing the nuclear reactor physics are based mainly on two theoretical areas: neutron transport theory and neutron diffusion theory, where it is necessary to remark that neutron diffusion theory is really a simplification of the neutron transport theory.
Numerical methods are used to solve the partial differential equations representing the nuclear reactor physics, and these methods are derived from discretization techniques. For numerical solutions in any scientific area, computational tools have been developed including software and hardware. In the past, the former computer processing was the sequential execution of computer commands, meaning to say that program tasks are carried out one after one. Modern computational tools have been developed for parallel processing, executing several tasks concurrently.
The computing branch dealing with the system architecture and appropriate software related to the simultaneous execution of computer instructions and applications is known as parallel computing science. Former developments in parallel computing were made in the late 1950s, following the construction of supercomputers throughout the 1960s and 1970s. Nowadays, clusters are the workhorse of scientific computing and are the dominant architecture in data centers.
Since the late 1950s, the performance of safety analyses was essential in the nuclear industry, in research reactors, but mainly safety analyses of nuclear power plants for commercial purposes. Scientific computing calculations were vital to these safety analyses, but with important limitations in computer/computing capabilities. At the beginning, the objective was to give a solution to partial differential equation models based on neutron diffusion or neutron transport with technology and methods available in those years. Numerical techniques were used first with finite differences and finite element approaches, and gradually up to now, with nodal finite element methods (NFEMs). Despite the numerical method employed, the computer code user faces the problem of solving extremely large algebraic systems challenging hardware/software capabilities. Generation of results for any reactor simulation in considerable short times is a desirable achievement for computer code users .
Recent developments of high-performance computer equipment and software have made the use of supercomputing in many scientific areas possible. The appropriate selection of parallel computing software, like newly developed linear algebra libraries, to be used in a specific project may result in a suitable platform to simulate nuclear reactor states with relatively prompt results.
Throughout the world, several research projects in the last decade have been developed with the main objective of making full tridimensional (3D) coupling simulations of nuclear reactor cores, leaving aside the obsolescence of the point kinetics theory. Most of the modern nuclear reactor simulators are based on neutron transport theory, or on neutron diffusion theory, to obtain detailed 3D results. As light water is used for cooling/moderating light water reactors (LWRs), a comprehensive analysis of the reactor core physics must include thermal-hydraulic phenomena, so that modern simulations perform reactor calculations with thermal-hydraulic feedback coupled with neutron kinetics calculations.
All the discussions included in this chapter are centered in a simulator for light water reactors. The computer code AZtlan KInetics in Neutron Diffusion (AZKIND) is part of the neutronic codes selected for their implementation in the AZTLAN Platform1 project in which neutron transport and neutron diffusion codes are being developed in Mexico. A (TH) model has been implemented recently and coupled with the neutronic (NK) model, and both models are based on HPC implementations.
2. Reactor core calculation overview
Although there has been growing interest in the transport-based core neutronics analysis methods for a more accurate calculation with high-performance computers, it is yet impractical to apply them in the real core design activities because their performance is not so practical on ordinary desktop or server computing machines. For this reason, most of the neutronics codes for reactor core calculations are still subject to the two-step calculation procedure, which consists of (1) homogenized group neutron parameters generation and (2) neutron diffusion core calculation.
In the core calculation steps that are the main concern of this work, nodal codes based on the diffusion theory have been used to determine the neutron multiplication factor and the corresponding core neutron flux (or power) distribution. Practically, almost all nuclear reactor simulation codes employ the two-group approach involving only fast and thermal neutron energy groups for the applications to light water reactors (LWR). However, numerical calculations with the two-group structure are not appropriate in the analysis of cores loaded with mixed oxide fuels or analysis of fast breeder reactors, since the neutron spectrum is influenced more by the core environment, requiring much more energy groups than only two groups.
As settled in Ref. , even using a high-performance computer, a direct core calculation with several tens of thousands of fuel pins is difficult to perform in its heterogeneous geometry model form, using fine groups of a prepared reactor cross-section library. The Monte Carlo method can handle such a core calculation (see also the Serpent code), but it is not easy to obtain enough accuracy for a local calculation or small reactivity because of accompanying statistical errors, besides the large calculation times. Instead of using neutron transport computer codes, the nuclear design calculation is performed in two steps: (1) lattice calculation in a two-dimensional infinite arrangement of fuel rods or assemblies for the generation of homogenized lattices jointly with their corresponding homogenized cross-sections and (2) core calculation in a three-dimensional whole core, with a neutron diffusion code using the information of the previous step.
As shown in Figure 1 , the lattice calculation prepares few-group homogenized cross sections which maintain the energy dependence (neutron spectrum) of nuclear reactions, and these reduce the core calculation cost in terms of time and memory. The final core design parameters are not concerned with continuous energy dependence, but spatial dependence, such as power distribution, is important to avoid high local neutron fluxes or high absorbing materials causing significant neutron flux gradients, mainly when safety analyses are performed upon the final proposed core designs.
In the core calculations with space-dependent data (cross sections and neutron flux), the effective cross sections are processed, with a little degradation in the accuracy as possible, by using the results from the multi-group lattice calculation. Lattice code calculation and codes are not discussed here.
There are two processes followed for lattice calculation. One is the homogenization to lessen the space-dependent information and the other is group-collapsing to reduce the energy-dependent information as shown in Figure 2. The fundamental idea of both methods is to preserve neutron reaction rate. The next step is to consider the conservation of reaction rate in the energy group
The number of few groups depends on reactor type and computation code. Two or three groups are adopted for the NK- and TH-coupled core calculation of LWRs and much more groups (18, 33, etc.) are used for the core calculation of LMFRs (Liquid Metal Fast Reactors). Currently, revised methods exist for the improvement of cross-sections generation using computer codes dedicated to lattice calculation for few-groups approach, like in Ref. , where three topics are involved: (1) improved treatment of neutron-multiplying scattering reactions; (2) group constant generation in reflectors and other non-fissile regions, leading to the use of discontinuity factors in neutron diffusion codes; and (3) homogenization in leakage-corrected criticality spectrum, in which several leakage corrections are used to attain criticality, accounting for the non-physical infinite-lattice approximation. Another improvement was done in Monte Carlo codes , implementing reliable multi-group cross-sections calculations for collapsed flux spectrum. Ref.  focuses on calculating scattering cross sections, including the group-to-group scattering.
The following sections contain, as a matter of example, summarized explanations of the AZKIND nuclear reactor simulator in which the reactor physics is based on neutron diffusion theory.
3. Neutron diffusion theory and nodal methods
3.1. Multi-group time-dependent neutron diffusion equations
In addition to boundary conditions for neutron fluxes, initial conditions must be satisfied by neutron fluxes and neutron precursor functions. Parameters involved in the above equations are described in .
3.2. Spatial discretization
The spatial discretization of Eqs. (1) and (2) is strongly connected with the discretization of a nuclear reactor core of volume Ω. Representing the neutron flux and the precursor concentrations in terms of base functions defined over Ω, it is possible to write
where and Table 1 contains the expressions representing the calculation of each matrix coefficient.
3.3. NFE method in spatial discretization
As fully explained in  and summarized in , a simple NFE element is characterized by the fact that for each node, the function unknowns to be determined are the (00) Legendre moment (average) of the unknown function over each face of the node and the (000) Legendre moment over the node volume. Figure 3(a) shows a physical domain Ω graphically represented after generating an
In the NFE method RTN-0, the normalized zero-order Legendre polynomials defined over the unit cell Ω
The matrix elements are quantified introducing the following nodal basis functions :
An extensive discussion on nodal diffusion methods can be found in Ref.  for space discretization using simplification approaches for calculating the moments over a node.
3.4. Discretization of the time variable
Once the spatial discretization is done, the
where , and
For time integration, parameters
Once the formulation to be used for time integration is established, the
For a known vector the algebraic system (9) is solved for the neutron fluxes . Therefore, the computing process requires an initial flux vector for the first time step, which is used in (9) to determine new neutron fluxes at the end of the time step, thus using these neutron fluxes to calculate a new delayed neutron precursor concentration vector. This process is sequentially performed for each time step over the total time interval (
4. Reactor power distribution
Once the computer model to solve the reactor kinetics Eqs. (1) and (2) is able to provide the neutron flux profile, the next objective is to know the power distribution in the reactor configuration. It is necessary to be aware that the neutron flux is by itself the shape of the power distribution in multiplicative materials. The numerical methods presented in previous sections to solve Eq. (9) produce an algorithm capable to obtain the neutron flux profile for a reactor steady state. The calculated neutron flux has the following property over the domain Ω: . To determine the real average neutron flux in the reactor core, , it is necessary to specify the magnitude of the fluxes. For instance, a flux normalization factor can be introduced such that
Theoretically, it would be best to determine the flux level which resulted in a critical reactor . This could be accomplished by coupling of the NK model with the TH model of the whole reactor. In practice, however, the scaling factor is determined such that the total generated thermal power corresponds to some user-specified value
where is the volumetric heat generation rate in the fuel in units of [W/cm3],
where is a dimensionless factor, is the energy released by a nuclear fission reaction in [MeV/fission], and the sum over is the volumetric fission rate in [fissions/(cm3∙s)]. Thus, Eq. (10) is written as
In a more general way, for a reactor volume
Therefore, using the reference total thermal power specified by the code user, the flux normalization factor can be written as
where the factors “kappa-fission” are With the flux normalization factor calculated as above, the actual thermal power distributions in the reactor core can be calculated using the current neutron flux in the reactor core . Nevertheless, it is necessary to introduce the value of . This value is used as an average energy released of 200 MeV (i.e.,), based on the energies released by the fission of the U235 nuclei .
In summary, once the NK model is used to generate the neutron flux distribution in the reactor core, expression (12) can be used to calculate the thermal power being generated along all the nodes in a thermal-hydraulic channel of area and height H. This thermal power can be the axial power profile needed by the TH model to produce the thermal-hydraulic state corresponding to the generated thermal power.
5. Neutronic and thermal-hydraulic coupling model (NK-TH)
The description contained in this section is based on a work published by Ceceñas in Ref.  about a TH model developed for boiling water reactors. The TH model was modified from a point kinetics approach with an extension of the NK model to 3D and implemented in the development of AZKIND.
The treatment of neutron kinetics in  has been improved by coupling a 3D solution of the neutron diffusion equations with an arrangement of TH channels in parallel. Each channel independently contemplates three regions: (1) one phase, (2) subcooled boiling, and (3) bulk boiling. The objective was to implement a detailed model of a nuclear reactor core, which is somehow perturbed to simulate NK-TH coupling. These perturbations are obtained when the power generated in a group of channels changes and thus affecting the TH state of each channel.
The original  TH model is based on a generic channel, which is adapted by transferring to it the operational data as flow area, generated power, axial power profile, and subcooling, among other parameters. Each channel is associated with a number of nuclear fuel assemblies and an axial power profile. Although the neutron model is a two-dimensional model for the radial power profile in each
For the implementation in AZKIND of the TH model of Ceceñas, the grouping of fuel assemblies was maintained for generating a reduced number of TH channels; operational data are also used. The main difference is that the NK model recursively computes the axial power profile for each channel, and this thermal power is the updated source of power for TH model. Therefore, a “new” thermal-hydraulic condition is generated, and it is used by the NEMTAB model to update the nuclear data to generate new thermal power profiles with the NK model. The process is iterative, and it stops when the convergence is met. Convergence is achieved when updated conditions do not change in both NK and TH models.
The NK-TH coupling in AZKIND performs core calculations as described above to obtain a steady-state reactor core condition. For transient conditions in a time interval
The TH model comprises the solution of the mass, momentum, and energy conservation equations in the three regions contemplated by the channel: (1) one phase, (2) subcooled boiling, and (3) bulk boiling. The system receives heat through a non-uniform source whose profile is axially defined plane by plane. This axial use of the power profile allows the inclusion of a wide range of axial profiles, from relatively flat to profiles with their peak value at some axial point in each channel in the core.
5.1. Heat transfer in the fuel
The heat transfer and temperature distribution in the fuel and cladding can be calculated by a simple model where the heat diffusion equation is solved in one dimension (radial) for a fuel rod, since the conduction in axial direction is small compared to the radial one, it can be neglected. An energy balance per unit length yields
where and represent thermal resistances per unit length. The coefficient of heat transfer to the refrigerant fluid is calculated by the Dittus-Boelter or Chen correlation, depending on the type of flow, which can be in one or two phases. These equations are used for the radial averaging of the temperatures in the fuel rod.
5.2. Reactor coolant dynamics
The conservation equations of mass, energy, and momentum are applied in this case to a flow of water along a vertical channel, where the dynamics of the fluid heated by the wall of the fuel is modeled. Conservation equations can be expressed as 
In this work, the conservation equations are solved by the Integral Moment method , according to which it is assumed that the refrigerant is incompressible but thermally expandable, and the density is a function of enthalpy at a constant pressure
Neglecting terms related to pressure changes and wall friction forces, the energy equation is simplified as
where the axial flow variation can be obtained by
This equation provides the flow variations with respect to an average value imposed as a boundary value or provided by the dynamics of the coolant recirculation system. Three regions are defined by which the coolant circulates as it ascends into the channel: a one-phase region, a subcooled boiling region, and a bulk boiling region. The first region begins at the bottom of the channel, where the coolant enters with known enthalpy and ends at the point of separation of the bubbles
where, in this case, the parameters ks and
ks = 0.71 + 1.2865 × 10−3
The total pressure drop in the channel is made up of the contributions of each region. Every term in each region includes the contribution by acceleration, gravity, and friction. For the channel arrangement, the steady state is obtained by iterating the coolant flow rate of each channel to obtain the same pressure drop for all of them. This iteration consists of a correction to the flow defined by the deviation of the pressure drop of the channel with respect to the average of all the channels:
where is the flow rate for channel
It is observed that even though the pressures are equaled, the value of the pressure drop in the core is not imposed as a boundary condition. Convergence is achieved when the following relationship is met: . By changing the flow rate of the channel for each iteration, the enthalpy and void fraction profiles are affected. It is necessary to recalculate the TH solution at each iteration for all channels, achieving convergence when every parameter involved in the thermal-hydraulic calculation remains unchanged.
5.3. Neutron kinetics: thermal-hydraulics (NK-TH) coupling model
Although reference  has important issues to be considered in the development of an NK-TH-coupled model, those issues are not repeated here, but taken into account. The most direct way of coupling NK module and TH module, as implemented in AZKIND, consists simply in that axially both NK mesh and TH mesh have the same partition, making possible to assign an NK node at position
As it can be seen in Figure 4, before initiating the NK-TH feedback process, the initial nuclear parameters and kinetics parameter (XS) are loaded from files constructed in NEMTAB format, previously generated by means of a lattice code. Then, following the reading of the nuclear reactor burn-up state and thermo-physics initial conditions, the XS parameters are obtained from the Nemtab multi-dimensional tables by means of interpolation calculations.
The process continues as follows. The corresponding neutron flux is calculated in the NK module with the
The important variables sent to the NK module are the fuel temperature (
6. High-performance computing in AZKIND
6.1. PARALUTION linear algebra library
HPC was implemented in AZKIND with the support of the linear algebra solvers library PARALUTION . This open-source library is optimized for parallel computing process using graphics processing units (GPUs). For the numerical solution of an algebraic system PARALUTION includes numerical solvers to obtain the solution vector for a known vector and a specific matrix that can be a symmetric or a non-symmetric matrix being also a sparse or a dense matrix. The working matrices in AZKIND are sparse non-symmetric matrices, and the
6.2. Parallel processing for neutronic model
To demonstrate the HPC implementation in AZKIND, as described in Ref. , very large matrices were constructed for fine spatial discretization of arrangements of nuclear fuel assemblies of an LWR. Fine discretization means that each fuel assembly was subdivided in a mesh of size 10 × 10. As an example, an arrangement of 6 × 6 fuel assemblies consists of a square with 36 fine-discretized fuel assemblies. The corresponding algebraic system for each fuel arrangement was solved with parallel processing performed by the
Matrix dimension (
|1 × 1|
|2 × 2|
|4 × 4|
|6 × 6|
|10 × 10|
|GTX 860 M||2.1||7.9||31.3||No memory||No memory|
|Tesla K20c||1.3||4.0||16.6||40.1||No memory|
|GTX TITAN X||1.0||2.6||10.4||26.7||95.4|
|Assemblies||1 × 1||2 × 2||4 × 4||6 × 6||10 × 10|
|GTX 860 M||11||16||12||–||–|
|GTX TITAN X||24||48||36||37||26|
Figure 5  shows the distribution of nuclear fuel assemblies in the core of a boiling water reactor. Excepting the blue-shaded zone, colors are for different types of fuel assemblies. In the plane
As described in , a reactor power transient was simulated as the capability to remove neutrons was highly increased in the perturbed assembly shown in Figure 5. An increase as step function in the neutrons removal capability during 3 s is implemented in the perturbed assembly, after that the perturbation finishes and the transient lasts for two more seconds, giving a reactor power reduction. The time step used in this simulation was 0.1 s. Figure 6 shows the power behavior over time, departing from a normalized value of 1.0 and reducing the power reactor to almost 80% of its original value. This reactor power transient was simulated with the AZKIND code, running on the three different GPUs listed in Tables 2 and 3. The right side of Figure 6 shows the time spent by AZKIND in a logarithmic scale, running in a sequential mode (Serial bar) and the times spent by each GPU card.
7. Simulation of a reactor core condition
A simple example was prepared to show the capability of the AZKIND code running with NK-TH coupling, and the thermal-hydraulic effect on power distribution is compared to the power distribution resulted from the NK model running standalone.
This example was prepared for a two energy group, that is, fast neutrons and thermal neutrons. In LWR, the nuclear fissions of the fuel atoms are mainly coming from the thermal neutrons present in the reactor core. The effect observed in Figure 7 is that the TH feedback induces an increase in the thermal neutrons population and so increasing power. As the coolant/moderator enters the reactor core through the bottom part of the reactor and the core is beginning the production cycle, the core design allows more power generation in the first third of the core active fuel. Also, as it was expected, in the map of fuel assemblies of the reactor core, the location of the fuel assembly with the highest generation of thermal power remained unchanged with the insertion of TH feedback.
8. Some advances on nuclear reactor simulation
In the last two decades, there have been significant advances in the development of nuclear reactor codes for 3D simulation with coupling NK-TH, supported with new modeling techniques and modern computing capabilities in software and hardware. Some examples of these advances are listed subsequently:
DYNSUB: Pin-based coupling of the simplified transport (SP3) version of DYN3D with the sub-channel code SUBCHANFLOW. See [16, 17]. The new coupled code system allows for a more realistic description of the core behavior under steady state and transient conditions. DYNSUB has successfully been applied to analyze the behavior of one eight of a PWR core during an REA transient by a pin-by-pin simulation consisting of a huge number of nodes. Some insights are pointed out on the convergence process with a detailed coupling solution modeling neighbor sub-channels and modeling adjacent assembly channels.
DYN3D: The code comprises various 3D neutron kinetics solvers, a thermal-hydraulics reactor core model, and a thermo-mechanical fuel rod model, see . The following topics are delineated in the reference: the latest developments of models and methods, a status of verification and validation; code applications for selected safety analyses; multi-physics code couplings to thermal-hydraulic system codes, CFD, and sub-channel codes as well as to the fuel performance code TRANSURANUS.
TRACE/PARCS: See . The study of the coupling capability of the TRACE and PARCS codes by analyzing the “Main Steam Line Break (MSLB) benchmark problem,” consisting of a double-ended MSLB accident assumed to occur in the Babcock and Wilcox Three Mile Island Unit 1. The model TRACE/PARCS generated data showing that these codes have the capability to predict expected phenomena typical of this transient and the related NK-TH feedback.
COBAYA3: See . This reference describes a multi-physics system of codes including the 3D multi-group neutron diffusion codes, ANDES and COBAYA3-PBP, coupled with the sub-channel thermal-hydraulic codes COBRA-TF, COBRA-IIIc, and SUBCHANFLOW, for the simulation of LWR core transients. Implementation of the PARALUTION library to solve sparse systems of linear equations was done. It features several types of iterative solvers and preconditioners which can run on both multi-core CPUs and GPU devices without any modification from the interface point of view. By exploring this technology, namely the implementation of the PARALUTION library in COBAYA3, the code can decrease the solution time of the sparse linear systems by a factor of 5.15 on GPU and 2.56 on a multi-core CPU using standard hardware.
CNFR: See . This reference summarizes three methods, implemented for multi-core CPU and GPU, to evaluate fuel burn-up in a pressurized light water nuclear reactor (PWR) using the solutions of a large system of coupled ordinary differential equations. The reactor physics simulation of a PWR with burn-up calculations spends long execution times, so that performance improvement using GPU can imply in a better core design and thus extended fuel life cycle. The results with parallel computing exhibit speed improvement exceeding 200 times over the sequential solver, within 1% accuracy.
9. Conclusions and remarks
The state of the art in the topic of nuclear reactor simulations shows significant advances in the development of computer codes. A wide range of applications focusing, besides on improving nuclear safety, on more efficient analyses to improve fuel cycles/depletion have been found in a recent study. A considerable “saving time” factor in obtaining nuclear reactor analyses has been observed.
One important part of a nuclear reactor simulator is the benchmarking process to demonstrate reliability and repeatability in the simulation of real cases, for which data from reactor operation or comprehensive data from experiments are well documented. In this sense, extensive documentation is necessary for theoretical basis, numerical techniques and tools, and validation of both codes and simulation models.
- This work was performed under the auspices of the financial support from the National Strategic Project No. 212602 (AZTLAN Platform) as part of the Sectorial Fund for Energetic Sustainability CONACYT—SENER, Mexico.