A TOPOLOGY OPTIMIZATION PROCEDURE FOR ASSISTING THE DESIGN OF NUCLEAR COMPONENTS

In this paper, we show that a module implemented in the MCNP transport code to perform sensitivity analyses can be diverted to perform topology optimizations of nuclear equipment. Component design with this approach leads to sophisticated solutions that outperform their human-designed counterparts.


INTRODUCTION
Suppose the operators of a silicon doping unit wish to design a neutron moderator that homogenizes the distribution of neutron capture sites in an irradiated ingot.Or suppose a medical physics team wish to design a gamma-ray collimator, which has to be as efficient as possible in a limited volume.At present, these design studies are of parametric nature.The characteristics of the aforementioned components are optimized using an intuitive subset of parameters, which includes the choice of materials to be used, as well as some geometrical parameters, e.g. the dimensions of subparts, the aperture angles, and so on.For each combination of these parameters, a simulation of the candidate device is performed, and the best solution emerging among them is chosen according to the objectives of the designers.This design practice is standard but it has a weak point.It relies on human intuition, which provides with a general idea of the solution but generally lacks precision.To outmatch this human-based design, a transdisciplinary approach called topology optimization is blooming, alongside with the increase in computing power [1].This discipline aims at automatizing the design of equipment using algorithms able to take into account the physical laws involved while managing the objectives and constraints of the design.These algorithms already yield subtler, more performant solutions than human-based designs in a number of industrial fields [2][3][4][5].In our field, however, there are still seldom applications of these techniques, presumably because of the large computing cost associated with an accurate solution of the Boltzmann equation.At the present times, the most advanced studies are focused on the development of radiation shields for spatial applications and on the fuel organization problem [6][7][8].In these proceedings, relying on the theoretical framework proposed in [9], we will show that the PERT module implemented in the MCNP code [10][11][12] is precise and fast enough to perform realistic topology optimization of nuclear components.We will illustrate this progress with an application of practical interest, the design of a neutron spectrum shaper, and propose a benchmark on the optimization of the fuel distribution in a critical assembly.

SOLVING A TOPOLOGY OPTIMIZATION PROBLEM WITH THE MCNP CODE
Consider a block of matter, whose mass density at a point r = (x, y, z) of space is denoted (r).This block can be machined in various materials, whose respective atomic concentrations at the point r form a vector, denoted (r).A source, which can be external or intrinsic, fires particles into this block, whose transport can be modelled using a linear Boltzmann equation, ,(r,E,,,) = Q(r,E,), provided the density of particles remains low enough so that they do not interact with each other.In this equation,  denotes the angular particle flux, function of the particle energies and directions, E and ,  is the Boltzmann operator, which depends on the physical characteristics  and  of the block, and Q is the distribution in space, energies and directions of the source particles.
For industrial or scientific applications, one could search the best physical properties, (opt, opt), of the aforementioned block of matter so that the particle population in a space region V optimizes a property desired by the operators of the particle source.For example, one could wonder how to machine a block of polyethylene so that the total flux of particles in V be as small as possible (design of a particle shield).To answer this kind of problem, one must solve an optimization problem, which can be written In (1), O is a functional acting on the angular particle flux .The command min/max O is the objective desired by the operators of the source: it could be e.g. the minimization of the total particle flux at a point r d for a radiation protection problem, in which case one would take min O = ∫∫(r d ,E,)dEd.It could be the maximization of the radiation dose delivered in an area V of a body for a cancer treatment problem, in which case one would take max O = V -1 ∫∫∫D(E)(r,E,)drdEd for r ∈ V, with D(E) being the flux-to-dose conversion function and V the volume of area V.The optimization problem min/max O is subject to some constraints, denoted C(,,(,)) = 0, the first one being the requirement the particle transport obeys the Boltzmann equation, and the other ones being additional constraints, for example a limitation over the maximum quantity of matter available (weight constraint), or a constraint over the k eff of the system.This formalism is further discussed in [9] section 1.To solve the optimization problem (1), one can start from its Lagrangian formulation where  is a vector of Lagrange multipliers.The optimal properties of the particle propagator, (opt, opt), which solve problem (1) are then given by the following system of equations, To solve this system, a solution is to segment the block into a large number N of small voxels,  i=1…N .In each voxel i, the density (r) and fractions (r) are assumed constant, equal to i and i.The derivatives (3) over functions  and  thus become derivatives over parameters i and i, whose computation can be performed with the PERT module implemented in the MCNP transport code [13].Practical instructions on how to use the PERT module for computing derivatives are given in [9] section 1 for the computation of ∂L/∂i and in [9] section P (within the supplementary material file) for the computation of ∂L/∂i.The PERT module relies on the differential operator sampling method, which is faster and more precise than the usual procedure based on the adjoint problem [1,7].With this tool, all derivatives (3) can be computed in a single, direct MCNP calculation, which makes the precise resolution of a type (1) problem feasible in a humanly compatible time with the computer resources available in the 2010s.Equations (3) are then solved using an iterative algorithm, described in [9] section 2.2.This algorithm starts from an uniform configuration, i = 0 and i = 0 ∀i, then improves, iteration after iteration, the design of the structure by incrementing the parameters i and i by small, discrete quantities, ± and ±, according to the constraints of problem (1).This resolution procedure has been successfully tested over a large range of problems [9].In these proceedings, we will illustrate its efficiency by solving two problems: the design of a neutron spectrometer in section 3, and the minimization of a critical mass in section 4.

EXAMPLE OF APPLICATION: DESIGN OF A NEUTRON SPECTRUM SHAPER
Of all the possible applications of the optimization procedure described section 2, the following problem stands out for two reasons: (i) it has a complex, non-linear objective functional O, whose handling is nonetheless possible with the approach described above; (ii) its formalism is transposable to a large category of problems, in fact all of those involving the calculation of a distance to an objective (see e.g. the design of a gamma-ray collimator in [9] section 5.2).This problem, therefore chosen to illustrate the methodology described in this paper, is stated below.
Problem.A 14 MeV neutron source is positioned at a point rs, at a distance 100 cm from the center of a detection area (DA).The operators of the source have in stock 20 metric tons (∼1.76 m 3 ) of natural lead, and would like to shape it so that the energy spectrum (E) of neutrons in the detection area be as close as possible to an objective spectrum obj(E).Find the design of such a neutron spectrum shaper.
First, let's tile space with N voxels i, and denote  the vector of the lead mass densities i=1…N in i, and V the vector of the volumes Vi of these voxels.The detection area (here a cylinder of length 12.5 cm and radius 6.25 cm) is surrounded by a layer 1 cm thick of natural boron, whose mass density b is also to be optimized.With these notations, the design problem can be formulated as follows In (4), P is the weight of the device, Pmax the maximum permissible weight (20 tons here), max the natural mass density of lead (11.34 g/cm 3 ), and d(,obj) the distance between the energy spectrum (E) in the DA and the objective spectrum obj(E), with w(E) being a weight function (here taken equal to 1/obj 2 ).The distance d is 0 iff (E) ∝  obj (E), which is the objective of the design.In [9], the author attempted to solve this problem without the weight constraint, and with a small number of voxels.The result showed promises, but was not fully satisfactory as questions arose on the usefulness of the outer regions of the calculated structure.The weight constraint is introduced here to rationalize the positioning of the material in the device.For this example, the chosen objective spectrum is typical of a fast reactor spectrum [14], and is displayed in gray on the right side of Fig. 1.On the left side of Fig. 1, the densities  in the voxels  obtained with the optimization procedure are indicated for several iterations of the algorithm, using a gray scale whose unit is g/cm 3 .The units of the y and z axes are cm.The full 3D structure of the neutron spectrum shaper is generated by rotating these density maps around the symmetry axis y = 0.The small red cell on the left houses the 14 MeV source, and the larger red cell on the right is the DA.For each map, the normalized spectrum obtained in the detection area is given on the right.The computations were performed using MCNP with cross-sections taken from JEFF-3.1.They lasted ∼3 days/iteration on 24 CPU for 10 9 source neutrons and 1145 cylindrical voxels.The evolutions of the distance d(, obj ) and of the weight P of the device with the iteration number n are given in Fig. 2. As can be seen in Fig. 1, the algorithm begins by accreting matter along the axis y = 0, which connects the source to the detection area, in order to slow-down the 14 MeV neutrons.A reflector is simultaneously assembled around this central hub to increase the number of neutron collisions in the structure and shift the spectrum towards lower energies.The resulting device is efficient: the objective spectrum is correctly imitated over 6 orders of magnitude in energy, starting from n = 35.Differences still remain between the final spectrum at n = 35 and the objective spectrum, especially at high and low energy.These differences may be reduced by using a composite material, made of a mixture of lead and lighter isotopes, which will allow access to different neutron slowing-down modes.

DESIGN SENSITIVITY TO THE TRANSPORT DATA
The resolution of problem (1) with a transport code naturally raises the question of the sensitivity of the computed design to the uncertainties on the transport data.More particularly, one may wonder if a slight modification of a transport datum, e.g. a change in the position or width of a cross-section resonance, could strongly alter the obtained optimal characteristics, ( opt ,  opt ), of the device.In other words, could problem (1) be subject to topological phase transitions?Despite its major interest, answering this question should prove difficult.Indeed, the resolution of the optimization problem is iterative, yet predicting the outcome of a complex iterative procedure is an open question in mathematics.To answer it, for a moment it will be necessary to be content with a Full Monte-Carlo approach.In this approach, a number of transport data sets could be sampled using the uncertainty data provided in the transport libraries, then for each data set problem (1) could be solved.In the end, the uncertainty on the optimal design parameters ( opt ,  opt ) could be evaluated by computing their average and variance.This approach is straightforward but, as the resolution of a type (1) problem typically requires times of the order of several days to several months.CPU depending on the number of voxels in the design, conducting it would require years.CPU, which is a computing power the author did not have access to at the time of writing these proceedings.In this section, a simple benchmark problem is hence proposed, which can make a good starting point for evaluating the uncertainties on the design arising from the uncertainties on the transport data, but also for comparing the designs obtained with other transport codes than MCNP and other optimization approaches than the one introduced in section 2. This problem is stated below.
Problem.Consider a spherical assembly, 50 cm in radius, containing a mixture of low-enriched uranium (LEU, chemical formula: 97% at. 238U + 3% at. 235U) and polyethylene (PE, chemical formula CH2).The assembly is made of 100 concentric spherical crowns  i , comprised between the radii R i and R i+1 , with R i = 0.5i cm.Each crown  i contains a fraction  i of LEU and a complementary fraction (1− i ) of PE.With these notations, the LEU mass M in the assembly is given by where  LEU = 19.0445g/cm 3 is the LEU mass density and V the vector of the volumes V i of cells  i .The PE mass density is 0.94 g/cm 3 .Find the LEU fraction profile  opt that allows the assembly to achieve criticality with the smallest mass M of fuel possible.
This problem can be solved using the approach described in section 2. The full details of the procedure to follow are given in [9] section P (supplementary material).The calculations were performed using MCNP with an iterative fission neutron distribution (KCODE mode, with 50 passive cycles and 100 active cycles of 2.0 10 5 neutrons), with cross-sections taken from two libraries, ENDF/B-VII.0 and JEFF3.1 (1) .They lasted ∼20 min/iteration on 24 CPU.The LEU fraction profile that minimizes the critical mass is given in Fig. 3, as a function of the distance r from the center of the assembly, for these two databases.The evolutions of the fuel mass M and the reactivity rho = (k eff −1)/k eff of the assembly, expressed in pcm unit, with the iteration number n of the optimization algorithm are given in Fig. 4 for these two databases.The statistical errors on the reactivity values are ∼9 pcm regardless of n.In Fig. 3, it is observed that the topology optimization algorithm generates a two-zone device: a fuel zone, comprised between r = 0 and 44 cm, surrounded by a reflector made of pure PE, for r > 44 cm.The fuel zone houses a small void at its center (2) , from which the LEU fraction swiftly rises to reach a maximum of 3.4%, then decreases when r increases, to finally fall at 0 at r = 44 cm.Overall, the LEU fractions are low everywhere.By thus diluting the fuel in the PE, the optimization algorithm generates a structure that efficiently thermalizes neutrons: 96% of the fissions in the device are due to thermal neutrons (E < 0.625 eV).The resulting structure is very efficient: 105.5 kg of 3% enriched uranium suffice to reach criticality.No significant differences are observed, whether in the design or the value of the minimum critical mass, between the results obtained with ENDF/B-VII.0 and JEFF3.1 libraries, which constitutes a preliminary element of answer on the robustness of the design.In Fig. 3 is observed the occurrence of small heterogeneities in the LEU fraction profiles, which are not mere statistical fluctuations.Indeed, the LEU fractions were purposely changed by very small increments,  = ± 5.10 -3 %, between two iterations of the algorithm (hence the artificially large number of iterations to converge).This increment is small compared to the amplitude of the observed heterogeneities, ∼0.2%.Hence, one might wonder if these heterogeneities play a role in reducing the critical mass.To investigate this question, the ENDF/B-VII.0LEU fraction curve was smoothed using a Savitzky-Golay approach [15] while preserving the mass M. The reactivity obtained with the smoothed profile remained compatible with 0. The observed heterogeneities are therefore probably artifacts.The author sees three possible reasons for their occurrence: (1) a segregation induced by decorrelations of fission sites that may occur in a metersize object; (2) a consequence of the non-zero statistical errors on the k eff perturbations, k i = (∂k eff /∂ i ), computed with MCNP.Fig. 5 illustrates the evolution with the distance r of this error, given in percent, computed for the optimal LEU profiles found with the ENDF/B-VII.0 and JEFF-3.1 libraries.However, complementary calculations involving higher numbers of source neutrons and KCODE active cycles were performed, yet showed similar heterogeneities, indicating that these two first leads are likely not the full answer; (3) a numerical issue induced by the printing format of the k eff perturbations in MCNP outputs.Indeed, while all other perturbations computed by MCNP are printed in scientific format, the k eff ones stand out as being printed as integers in pcm unit.Therefore, when a keff perturbation is very small, lower than 1 pcm, yet non-null, it is actually printed 0 in MCNP output files, thus inducing a hidden numerical error.These technical considerations highlight the interest of benchmarking the results obtained in this section with other optimization approaches, e.g. with procedures that could use evolutionary algorithms or adjoint methods, as in [6][7][8].

CONCLUSIONS
In this paper, we showed that the PERT module implemented in the MCNP transport code, usually used to perform sensitivity analyses, can be diverted to automate the design of nuclear components.Two practical optimization problems were then solved using this tool, the design of a neutron spectrum shaper and the minimization of a critical mass, illustrating the potential of the proposed methodology.

Figure 1 .
Figure 1.Iterative design of a neutron spectrum shaper with the optimization procedure.

Figure 2 .
Figure 2. Evolutions of the distance d and of the weight P of the device with the iteration number n.

Figure 3 .
Figure 3. Optimal LEU fraction profile as a function of the distance r from the assembly center.

Figure 4 .
Figure 4. Evolutions of the fuel mass M and the reactivity rho in pcm with the iteration number n.

Figure 5 .
Figure 5. Evolution with r of the statistical relative errors on the keff perturbation values at n = 500.