DiffMod – statistical 2D simulation model of neutron propagation and moderation

1. Introduction

Simulations of neutron propagation and moderation are important to understand and optimize target – moderator – reflector assemblies in neutron sources that provide neutron spectra different to the primary neutron spectrum emitted from the neutron production target. For most applications, as well for thermal neutron sources for neutron scattering or analytics, as for epithermal neutron sources for BNCT, industrial neutron imaging, or isotope production, the neutron spectrum has to be shifted to lower energies compared to the primary neutron spectrum.

Monte-Carlo methods have been established very well to describe nuclear particle interactions and transport since the Manhattan project in the 1940s [3]. Today, program packages like MCNP [2], FLUKA [1] or PHITS [7] are in use to describe neutron production by fission, spallation, or nuclear reactions as well as the neutron propagation and the interactions with any materials in their way. These programs follow the path of every single particle from the source (e.g. a proton beam) along all interactions with matter, e.g. conversion of a proton into a daughter neutron, energy loss and directional change by elastic scattering with matter, or absorption of a neutron with associated emission of a γ photon.

The interactions of the particles with matter are described by cross-sections that are available in public databases, e.g. the ENDF database published by the IAEA [5]. On the basis of the energy of the particle impinging, the material at the place where the particle is located, and the cross-section values, the probability of a scattering or nuclear reaction event is evaluated. On the basis of a random number generator and the probabilities mentioned before, the program decides if the reaction takes place and either follows the reaction products or continues to propagate the original particle to the next location, where again interaction with some material is possible. Every single particle is tracked from the source through all interactions until it (and all daughter particles) have either been absorbed or have left the volume of interest.

In the limit of an infinite number of particle tracks, this approach yields a very realistic and reliable view of the particle field inside matter, but typically a high number of particles ( 10 8 or more) needs to be tracked to get a suitable description of complex phenomena, as e.g. neutron moderation inside of and extraction from a moderator material of limited size. In case a high resolution in particle energy, time evolution or in space is required, more particle tracks need to be computed to receive sufficient statistical quality inside the fine bins. Such a simulation typically takes several days of computing time on a workstation cluster, so that a parameter variation of some geometrical details or a systematic variation of moderator material composition is tedious, but the simulation results then have a quality that is comparable to experimental outcomes as long as the data bases used contain the correct cross sections.

As the particles are computed one after the other, the evolution of e.g. a neutron cloud after a proton pulse in a pulsed neutron source can only be reconstructed when the time evolution of all particles has been recorded with sufficient time resolution and sufficient statistical quality. There is no generic information available on the evolution of the ensemble of particles.

For a quick variation of geometrical parameters and for the deeper understanding of the development of the global time-evolution of the neutron ensemble during the moderation process, a complementary approach is required. We intend to describe the evolution of the ensemble in time by describing the motion, the scattering processes of neutrons with the moderator or reflector material, and the absorption of neutrons. For this purpose, we have started to develop the program package DiffMod (from the main actions DIFFusion and MODeration). The main model details and first simulation results are described in this paper.

2. Model

To be able to describe the evolution of the ensemble of neutrons in time with a purely statistical approach, it is necessary to have a complete overview of the entire ensemble at the point of time currently under consideration. To be able to restrict the memory and computing time requirements, some important simplifications are necessary.

The neutron cloud is described by the position, the direction of motion and the kinetic energy of all neutrons. A reasonable binning of the population is necessary to group the ensemble into a finite number of units that can be treated individually.

The spatial resolution is simplified to a 2-dimensional square array of typically 1 cm edge length. A higher resolution can be chosen if necessary to describe fine geometrical details. As the mean free path of neutrons above thermal energies is typically well above 1 cm, this resolution is suitable to describe sufficiently well all thermal moderator and reflector material assemblies. An example of the spatial resolution doubled to 0.5 cm is shown in Fig. 4 and discussed there. In case of neutron absorbing materials, the mean free path for neutrons with energy above 1 keV is still above 1 cm. For strong absorbers of thermal neutrons, as e.g. pure Cd or B₄C, typically thinner layer thicknesses than 1 cm are used. In that case, the cell’s material needs to be a mixture of the absorber and the structural material around, so that the absorption probability within the cell is the average over the volume. The absorption probability for thermal neutrons is so high that the diffusion of thermal neutrons out of such a volume can be neglected anyway, so that the coarsening of the absorber geometry is not a problem.

The third dimension (perpendicular to the plane under investigation) is assumed to be homogeneous, i.e. a neutron propagating perpendicular to the 2-dimensional plane that is resolved will just propagate into the same material as it comes from. The simplification to 2D turns out to be reasonably precise to describe the behaviour of a neutron cloud in the central plane of a structure that is several cm thick. In Section 4, a comparison between a 3-dimensional MCNP simulation of a moderator sphere and a 2-dimensional DiffMod simulation of a moderator circle is shown.

The propagation directions are restricted to the 8 main directions on the compass, which are the directions towards the nearest and next-nearest cells in the 2D square lattice, so that the neutron transport towards these cells is simulated.

The energy resolution is logarithmic with two levels per decade. The highest energy depends on the initial neutron spectrum, the lowest energy is 31.6 meV for thermal neutrons. For neutrons from high energies down to thermal, one can assume the interaction with matter to be a purely atomic interaction. The binding energy of the atoms in the matrix as well as molecular excitations can be neglected.

The simplifications used here forbid the use of the approach to describe cold neutron sources, as the molecular interactions in the moderator material become important below thermal energies.

The program treats the evolution of the ensemble from one timestep to the next. The following interactions are considered:

Neutron absorption time constants depending on the neutron energy and the material are calculated at every place and the neutron intensity is reduced according to the length of the timestep.

Neutron scattering and moderation probabilities are described by the elastic scattering cross sections. These are treated locally, as well. It is assumed that the direction of the neutron after a scattering event is completely random, i.e. all the 8 directions will be occupied with the same probability after any scattering event.

The scattering probability is proportional to the elastic scattering cross section of the material and to the velocity of the neutron. Above thermal energies, the energy change of a central collision of a neutron with a nucleus with the mass A follows Eq. (1) with the initial energy E and the final energy E ′ . (1) ln ( E ) − ln ( E ′ ) = 2 ∗ ln ( ( A + 1 ) / ( A − 1 ) )

The energy transfer is lower, if the collision is not central. On average, the change per collision in the logarithm of the energy is (2) ln ( E ) − ln ( E ′ ) = ξ

According to the table published in Section 4.7.3 of Ref. [6], ξ = 0.49 ∗ 2 ∗ ln ( ( A + 1 ) / ( A − 1 ) ) within an accuracy of 2% for all materials except hydrogen. The values for the typical hydrogen or deuterium containing moderator materials are listed there.

From these values, the program calculates the number of scattering events necessary to reduce the neutron’s energy by a factor of 10 , which is the energy ratio between neighbouring bins. From the scattering cross section, one can then calculate the average time for a scattering event to take place and for a “moderation step” towards the next lower energy level. From these times known, the probability that one of these processes happens within the current timestep is derived, and the respective fraction of the neutron population is transferred to the next lower energy level.

3. Program details

The DiffMod program is written in python3 using the GR framework [4] for graphical output. It has been tested in a Linux environment as well as in the Windows 10 operating system.

For the primary neutrons produced by a 70 MeV proton beam impinging on a Ta target, a MCNP simulation result [9] has been converted to the DiffMod input spectrum and directionality. The simulations are sensitive to the initial neutron spectrum because the mean free path of neutrons above 100 keV kinetic energy strongly depends on the energy of the neutron, so that the shape of the neutron cloud is influenced by the initial energy spectrum.

Currently, the cross-sections of typical moderator and reflector materials as H₂O, D₂O, polyethylene (PE), Be, graphite, Fe, and Pb as well as B as absorber and the mixtures of these materials are integrated. The import of the cross-sections of all isotopes available in the ENDF database is in preparation (see Section 5).

The timesteps start with 0.1 ns and are increased during the run of the program. With increasing timesteps, the travel distance of high-energy neutrons may exceed the distance between neighbouring cells or the moderation probability towards the next lower energy may exceed 40%. In this case, inaccuracies will occur. The program checks the population of the different energy levels and will allow the increase of the timesteps, if the population of the energy levels that are treated properly is above a limit that can be defined in the input file, typically above 90%. With this criterion, these inaccuracies only occur when the respective energy levels are already depleted.

The neutron production is assumed to happen at the time t = 0 inside the target. The pulse shape of the neutron pulses is a convolution of the timescales of the events in the moderator-reflector assembly and the shape of the proton pulse. This convolution can easily be done analytically during the post-processing of the detector data, e.g. using the McTal software [8]. Including an elongated proton pulse in DiffMod would strongly increase the necessary calculation time because the initial neutrons have to be treated with highest time resolution during the entire proton pulse length.

An input file allows the definition of the geometry (i.e. area and cell size of the area under investigation) and the materials in every cell. In addition, detector volumes (2D) and detector planes (1D) can be defined. After every timestep, the energy-resolved neutron flux is written to an output file. The spectra recorded at the detector planes are directionally resolved, which makes them useful e.g. to measure the efficiency of the neutron extraction through an evacuated beamline. The output format of the detector data is inspired by the meshtally of MCNP, so that it can be read and treated with the same software environment that is available for MCNP outputs. We propose to use the McTal software that has been chosen for the plots presented in this paper [8].

A video sequence is written that shows the evolution of the neutron field at four different energy levels: 1 MeV, 1 keV, 1 eV and thermal (31.6 meV). This video allows one to understand in detail the propagation and moderation of the neutrons in a structured moderator-reflector assembly.

On a standard Windows 10 laptop a simulation of a moderator-reflector assembly with 30 × 50 cells (i.e. 300 × 500 cm² area) in the time range up to 1 ms after the proton pulse takes well less than one hour.

4. Benchmark against MCNP

OPEN IN VIEWER

Fig. 1.

Inner part of the MCNP geometry (left) and DiffMod geometry (right) for the PE moderator – Pb reflector assembly. The surface detector 1 and the volume detectors 2 and 3 in DiffMod compare to the right surface of the target volume 1 and the volumes 900 and 901 in the MCNP model.

For a benchmark between a 3-dimensional MCNP model and a 2-dimensional DiffMod analogon, we have chosen a simplified target – moderator – reflector geometry as shown in Fig. 1. The model consists of a Ta target with 100 cm² surface area that is hit by a 70 MeV proton beam. The target is embedded at 5 cm depth inside a moderator sphere of 20 cm diameter. The DiffMod geometry imitates the central plane of the 3-dimensional model, i.e. a circle with 20 cm diameter as moderator and a 10 cm wide Ta line as target. For the purpose of absolute calibration, it is assumed that the 1 cm thick slice receives 10% of the protons impinging. The moderator is surrounded by a large reflector volume that is well thicker than the active range. In the case of MCNP, a 30 cm thick spherical shell around the moderator sphere was chosen. To restrict the calculation time, 5 to 10 cm thickness of the Pb reflector and 15 to 20 cm thickness of the Be reflector were chosen in the case of DiffMod.

There is an area detector covering the right surface of the target (volume 1) in the MCNP model which is imitated by the detector line 1, and the detector volumes 900 and 901 from Fig. 1 (left) are imitated by the detector planes 2 and 3. In both programs, the detector values are normalized by the area or volume, respectively. In the case of DiffMod, it is assumed that the 2D-slice under investigation is 1 cm thick, so that the results can be compared directly.

Here, we present some comparison results for two moderator – reflector combinations: polyethylene (PE) moderator with Pb reflector and D₂O moderator with Be reflector.

OPEN IN VIEWER

Fig. 2.

Snapshot of the neutron distribution in a PE / Pb moderator – reflector assembly 240 ns after the proton pulse. The four quadrants show the spatial distribution of 1 MeV neutrons (top left), 1 keV neutrons (top right), 1 eV neutrons (bottom left) and thermal neutrons (bottom right).

For the case of a PE moderator and Pb reflector, Fig. 2 shows a snapshot of the simulation run 240 ns after the proton pulse. One can clearly see that the moderation speed is much higher in the central PE region compared to the outer Pb reflector. The top left quadrant of Fig. 2 shows that some neutrons in the Pb reflector still are in the MeV range, while much more of them can be seen in the keV range in the top right quadrant. In contrast, the PE moderator region is already depleted in both of these ranges. The neutrons inside the PE moderator already start to populate the eV energy range, which can be seen in the lower left quadrant of Fig. 2. The neutrons have not yet arrived in the thermal energy range.

This is well in agreement with the moderation times down to thermal energies of 7 μs in PE and 2720 μs in Pb given in the table in Section 4.7.3 of Ref. [6]. The model of the average change of the logarithm of the neutron energy being constant with every collision described there yields the relative change in neutron energy being proportional to the square of the time evolved. Together with the moderation times, the model predicts an average energy of neutrons in PE of about 30 eV and in Pb about 1 MeV after 240 ns. Of course, most of the neutrons that stay in the Pb volume have already experienced some interaction with PE, so that the energy is strongly reduced compared to that value.

OPEN IN VIEWER

Fig. 3.

Pulse shape at the target surface for thermal neutrons (left) and neutrons in the keV energy range (right) in a PE moderator – Pb reflector assembly. The MCNP results are shown in dashed lines, the DiffMod results in solid lines.

Figure 3 shows the time evolution of the neutron field at the target surface for thermal and keV neutrons in comparison between MCNP and DiffMod. The discrepancy between the peak intensities as well as the decay times in both ranges presented is in the order of 10%, which is very good if we consider the strongly simplified model. In the keV range it can be seen that the maximum in the DiffMod model arrives about 10 ns later than in the MCNP model. This small deviation is caused by the initial neutron spectrum in DiffMod, that does not take into account primary neutrons with energies below 1 MeV.

Figure 4 shows the thermal pulse shape and decay at the two detector sites at some cm distance from the target. In both cases, the pulse height and shape are reproduced well by the DiffMod simulations. The reduced decay at detector position 3 due to the neutrons that have been moderated slowly inside the Pb reflector and diffuse back into the PE moderator is properly reproduced. The initial peak at the first appearance of thermal neutrons at about 2 μs is due to the fact that the transport at early times is strongly dominated by the transport of the fast neutrons before the first scattering event. In the case of DiffMod the dilution of the primary neutrons with the distance from the target is not properly modelled due to the transport in only 8 fixed directions. At later times (>30 μs) the transport is dominated by the diffusion of thermal neutrons, then this initial effect is properly averaged out.

OPEN IN VIEWER

Fig. 4.

Pulse shape at the detector positions 2 (in black) and 3 (in blue) for thermal neutrons in linear (left) and logarithmic intensity scale (right) in a PE moderator – Pb reflector assembly. The MCNP results are shown in dashed lines, the DiffMod results in solid lines. The green curve is a variation of the DiffMod result with 5 mm cell size.

The green curve in Fig. 4 shows the result of a DiffMod simulation with higher spatial resolution, i.e. 0.5 cm cell size. In this case, the initial effects are reproduced identical to the previously reported case, but the decay is described better. This is due to the fact that the transport is always assumed to happen from the centre of one cell to the centre of the neighbouring cell. During decay due to the diffusion of neutrons out of the moderator region, this may overestimate the expansion of the thermal neutron cloud. The improved description of the decay behaviour is being paid with a high price of computing time. The simulation with 0.5 cm spatial resolution takes about 9 times longer than the simulation with 1 cm spatial resolution, a factor 4 is due to the higher number of cells and a factor of about 2.3 is due to the smaller timesteps to avoid neutron motion on distances larger than the next-nearest neighbour cell distances.

In the case of D₂O moderator and Be reflector, the scattering probabilities are much weaker than in the case regarded before. If we use the same geometry as shown in Fig. 1 (right), we observe a strong deviation between the DiffMod results and the MCNP simulations. The peak intensities in the DiffMod simulation are lower, the peak is reached earlier and the decay speed is faster. The reason for this deviation is that the mean free path of the neutrons in Be is larger compared to Pb, and that therefore the diffusion out of the model area is quite probable.

For this reason, the model had to be adapted by increasing the Be reflector thickness by 10 cm at every side. The results are shown in Figs 5 and 6. Here, the thermal peaks are less high, but broader in time compared to the PE / Pb assembly shown in Figs 3 and 4. This is due to the fact that the confinement of the neutron cloud in D₂O is weaker due to the smaller scattering cross section compared to PE, but the lifetime is longer due to the lower absorption cross section. The effect of reduced confinement can also be observed in the almost negligible difference between the detector positions 2 and 3 in Fig. 6.

The agreement between the MCNP data and the DiffMod simulations is again very good. This allows one to use DiffMod as a quick tool for evaluating the efficiency of (thermal) moderator – reflector assemblies that do not strongly depend on their 3-dimensional structure, but can be approximated well by a 2-dimensional projection.

OPEN IN VIEWER

Fig. 5.

Pulse shape at the target surface for thermal neutrons (left) and neutrons in the keV energy range (right) in a D₂O moderator – Be reflector assembly. The MCNP results are shown in dashed lines, the DiffMod results in solid lines.

OPEN IN VIEWER

Fig. 6.

Pulse shape at the detector positions 2 (in black) and 3 (in blue) for thermal neutrons in linear (left) and logarithmic intensity scale (right) in a D₂O moderator – Be reflector assembly. The MCNP results are shown in dashed lines, the DiffMod results in solid lines.

Currently we are working on a routine that allows one to import the cross sections for all isotopes listed in the ENDF database [5].

In addition to the volume and area detectors already realized we will offer 1- and 2-dimensional meshtallys. The extension to a full 3-dimensional model is partly foreseen in the code, but will not be realized soon.

The current version of the program is published online at the link https://go.fzj.de/diffmod. A mailing list for information about updates is offered by the corresponding author.