Abstract
Fast quasi-monoenergetic neutrons can be produced by accelerating charged deuterons on tritium solid targets. Benchmark experiments were performed in many laboratories with intense D-T neutron sources. The aim is to validate the computational models and nuclear data for fusion applications. The detailed information on the neutron source term is highly important for the benchmark analyses. At present, the MCNP family of codes cannot explicitly model the D-T reaction for Deuterons in the KeV energy range. The physics for the D-T neutron production was modelled at ENEA (Italy) in the SOURCE and SRCDX subroutines to compile with the MCNP source code. Some improvements to the original subroutines were introduced. The differential cross-sections for the D-T reaction from the ENDF/B-VII library were built into the code. The relativistic approach was implemented for neutron kinematics. The new D-T neutron source model was applied to the MCNP5 simulation of the tungsten integral experiment performed at the OKTAVIAN facility. The uncertainty associated with the realistic D-T reactions was separated from the total uncertainty of the source term. The outcome of the benchmark analysis was an improvement in the quality of the computational model of the experiment.
1. Introduction
The first device for the acceleration of charged particles has been constructed by Cockroft and Walton in 1932, the same year when Chadwick discovered the neutron. Since the 1950s tritium was available as a target material to produce ~14 MeV neutrons by the He reaction, which is still the first candidate for plasma fusion technology. Thus, these neutron generators have been widely used to study the interaction of the neutrons with the structural materials of fusion reactors. Other applications of the device include neutron activation analysis, low cross-section measurements, and neutron therapy. The broad resonance of the D-T cross-section at the deuteron energy of 109 KeV permits high yields of fast neutrons to be obtained with low-voltage deuteron accelerators [1]. Integral experiments for fusion technology have been collected in the International Database for Integral Shielding Experiments (SINBAD), which is coordinated by the OECD Nuclear Energy Agency (http://www.nea.fr/) and Radiation Safety Information Computational Center (http://www-rsicc.ornl.gov/).
A lot of effort is devoted to computational techniques for neutron transport. A real source is not isotropic or monoenergetic. Primary neutrons are generated by deuterons which have interacted in the target and lost part of their original energy. The experimentalists usually provide the source specifications, but these specifications need a detailed investigation to determine parameter uncertainties, whose different contributions are difficult to separate out in a real experiment. Special attention is needed to define the anisotropy of the source and the neutron yield as a function of its flight directions, as well as the spread in the emitted neutron energies.
The MCNP family of codes cannot explicitly model the deuteron transport and neutron emission from the D-T reaction for incident deuterons with energies in the KeV range. At ENEA (Italy), a D-T source model was developed, which includes the effect of the deuteron slowing down inside a titanium-tritium target and the D-T reaction kinematics. The model is implemented in MCNP subroutines (e.g., SOURCE, SRCDX), which define the neutron source when compiled with the MCNP source code. The subroutines are available in SINBAD database.
The present work is a case study to identify the essential information which is needed for the source definition in computational analyses. Section 2 illustrates the approach to the D-T reaction kinematics from first principles. These results allow a refinement of the original ENEA source subroutines (described in Section 3). The upgraded D-T source model is internally validated by reproducing the theoretical angle dependence of the neutron energy and yields and by cross-validating the cell flux with the point detector flux. In Section 3 a brief description of the deuteron transport model is also provided. The sample case (Section 4) concerns the 14 MeV neutron facility OKTAVIAN, which was built at Osaka University and has been operating since 1981. In the framework of the SINBAD project the measured neutron leakage spectra from a Tungsten sphere were analysed. In this work the source model from the previous analysis is compared to the MCNP5 calculation with the improved neutron source model. The outcome of the benchmark analysis represents an improvement in the quality of the computational model of the experiment.
2. Theoretical Framework of the D-T Reaction
The He reaction is properly modelled as a two-body reaction. Some nomenclature: in Figure 1 the subscript “” denotes the incident particle (deuteron), “” the target particle (tritium), “1” the emitted neutron and “2” the alpha particle. and are the momenta in the laboratory frame and the angles between the outgoing particle and the deuteron direction. The energy of the outgoing neutron is calculated according to the relativistic equation (1). is the relativistic energy of the deuterons, which are given a nominal kinetic energy of 280 KeV: The effect of the relativistic kinematics formulation is studied in comparison with the classical kinematics equations. The relativistic equations are implemented in a MATLAB function (RELKIN.M) for graphical purposes (Figure 2). The original source for the particle masses is the NIST reference on constants, units, and uncertainties web site (http://physics.nist.gov/), except for triton. This is computed subtracting the electron contribution to the tritium mass taken from the atomic mass tables of G. Audi and A. Wapstra (http://www-nds.iaea.org/amdc/web/masseval.html). It is worthwhile to point out that in the ENDF-6 formats manual [2] this datum is incorrect. The calculation of the -value for the D-T reaction is 17.589 MeV. The relativistic calculation agrees with the one provided by the DROSG2000 code (http://www-nds.iaea.org/ndspub/libraries/drosg2000/), which also uses the relativistic approach. The classical results are obtained from the kinematics formulas available in the ENDF-6 formats manual. The difference in energy when solving the two-body reaction problem with the classical equations is about 30 KeV.
The angular distributions of the outgoing neutrons in the centre-of-mass system are represented by Legendre polynomials. The Legendre coefficients up to the tenth order are available in the ENDF files (MAT = 2, MF = 50, MT = 6) at the nominal deuteron energy of 280 KeV. For the conversion into the angular distribution in the laboratory frame it is necessary to calculate the momentum of the deuteron , whose evaluation is trivial, and of the neutron and apply the Jacobian as follows: is a constant independent of the angle, so it is adjusted to normalise the distribution at zero degree. This normalisation is a common practice because the angular dependence of both the neutron energy and the neutron yield on the angle is comparatively small in the forward direction [3]. The effect of both the anisotropy in the CM system and the transformation into the LAB frame can be assessed in Figure 3. For the CM anisotropic reaction the data of DROSG2000 code are replaced by the Legendre expansion from the ENDF/B-VII library. The same coefficients are used in the RELKIN function. The calculations with the new equations coincide with the DROSG 2000 results for the angular distribution in the LAB system. It is noticeable that the anisotropy of the reaction with 280 KeV deuterons implies a difference of ~3% at backward angles in the laboratory angular distribution. As a matter of fact, at energies typical for deuteron accelerators the assumption of isotropic CM distribution is to be considered a rough approximation.
3. The Deuteron Beam Slowing Down
The deuteron particles hitting solid Titanium-Tritium targets get slowed down by interactions with the electrons and nuclei in the medium. The electronic stopping power, that is, the inelastic scattering with electrons, can be considered independent on nuclear interactions. On interaction with Ti or T nuclei, the deuteron deflects and loses part of its kinetic energy by elastic collisions. Eventually, the deuteron interacts with the tritium nucleus to produce one neutron. For low-energy deuterons, the only reaction branches to consider are elastic scattering and nuclear fusion. The energy and angular distributions of incident deuterons inside the target cause the primary energy spread of the neutron energy at a fixed detection angle.
For the simulation of a D-T neutron source the first choice was to search into the MCNP family of codes. The deuteron ion transport in MCNPX is performed by physics models designed for incident deuterons of 200 MeV and higher. They are able to provide a good response in the 10–100 MeV range, however, the results are poor down to the MeV range. Till now, there are no D and T nuclear data libraries. Therefore, the D transport within the solid target can be achieved with MCNPX, but results are to be interpreted with caution and should be benchmarked against experiments or otherwise known answers. The D energy cut-off for the D-T fusion reaction in MCNPX is 2 MeV, so neutrons cannot be generated at the low energies of our interest. Starting from MCNP4 the possibility was introduced to develop and compile ad hoc source subroutines. The SINBAD database makes available the subroutines to compile with MCNP5 source codes for the simulation of the neutron emission in the Frascati neutron generator. The routines were provided by Angelone et al. from ENEA, who developed a detailed model of the FNG neutron source [4]. For the D slowing down inside the solid FNG target the basic-language source code from TRIM software was translated into FORTRAN language. The Monte Carlo program, the Transport of Ions in Matter (TRIM), is the most comprehensive program included in SRIM (http://www.srim.org/). SRIM is a well-established collection of software packages which calculate many features of the transport of ions in matter, such as ion stopping powers and ranges in targets, ion implantation, sputtering, ion transmission, ion beam therapy. The main author is James F. Ziegler. TRIM does not model nuclear reactions. M. Pillon implemented the D-T neutron production in the TRIM-based model for the D transport inside the Ti-T target. The essential characteristics of the ENEA D-T source model are the following:
(1)Monte Carlo method,(2)free flight paths between collisions, which “condense” negligible amounts of energy transfer and deflection angles,(3)impulse approximation for the free flight paths in the low-energy range,(4)universal interatomic potential,(5)Rutherford scattering at higher ion energies,(6)electron stopping power data from the TRIM code tables,(7)D-T double differential cross-sections retrieved from the DROSG2000 code calculations at different deuteron energies,(8)neutron generation by the modified von Neumann rejection method,(9)classical kinematics. The MCNP5 source subroutines are SOURCE.F90 (coding the D slowing down in solid Ti-T target and modelling the fusion reactions), SRCDX.F90 (specific for transporting neutrons generated in SOURCE to point detectors), and six other subroutines for numeric calculus. The D-T source routines require the use of the RDUM card in the MCNP input file to specify deuteron beam energy, target thickness, T/Ti atomic fraction, beam width, and target axis coordinates.
At ENEA the last version of the D-T source model was prepared for MCNP5 on the Unix systems. The possibility to perform routine calculations on Windows systems is achieved by producing a patch file to apply to the original MCNP5 source code for Windows systems.
Some flaws have been found in the last version of the source routines. The range of the deuterons was underestimated with reference to the TRIM calculations. The problem has been found in the stopping power of the Ti-T mixture. Since the units of the data imported into the routines from TRIM tables changed from the original MCNP4B version, some parameters needed to be modified accordingly. It is now possible to ascertain that the mean average range of the deuteron ions is about 1.5 μm, about the same value given by the TRIM code.
A new condition is implemented into the source routines to terminate the Monte Carlo history when the cumulative deuteron free flight path exceeds the target thickness.
The calculation of the neutron spectra on a very fine energy binning with the original source routines showed an unphysical spike at the deuteron maximum energy, corresponding to uncollided deuterons. After discussion with M. Pillon the special treatment of the first collision was dropped and the spike disappeared.
The major new contribution to the D-T source routines involves items 7 and 9 of the list above. The use of double differential cross-sections in the ENEA code had two drawbacks. First, it was not easy to check if the data were wrong or approximate because the data section was huge. Second, the cross-section data relied on DROSG2000 Legendre coefficients, which could change according to new experimental evidence or refinement of theoretical nuclear models. The source routines are modified to reconstruct the double-differential cross-sections from tabulated reaction cross-sections and the Legendre coefficiens, defined by DATA statements in the code with values as given in the ENDF/B-VII library. The results of a test calculation comparing the original source routines and the modified routines with the data from ENDF/B-VII are shown in Figure 4. The deuteron energy is 280 KeV. The other parameters are set to their standard values: the beam width is 0.7 cm, T/Ti atomic ratio is 1.5, and target thickness is 0.001 cm. The energy binning is 10 KeV. The source and the detectors are placed in the void. In the MCNP5 input file point detectors are located at different angles from 0 to 180 degrees. The slight differences between the cross-sections in last ENEA version of the source routines (labelled “enea”) and the modified routines with ENDF/B-VII cross-sections (labelled “leg.”) result in a small difference in the absolute yield (~5%). Except for this detail, the distributions at each angle agree very well.
For the present work the source routines are further improved to introduce the neutron relativistic kinematics as given by equations defined in the previous section. The use of relativistic kinematics results in a shift in energy, as seen from Figure 5.
The modifications above imply a revision of the SRCDX subroutine for point detectors. A cross-validation of the source model coding is performed by comparing the MCNP5 cell flux estimator (“cf”) and the point detector (“pd”) flux. Figure 6 shows that both tallies are consistent, except from some fluctuations in the cell flux estimator. The energy range starts at the theoretical neutron energy. The energy spread is greater in the forward direction. The shape of the neutron spectra can be explained by considering that the maximum of the He integral cross-section is at 109 KeV.
To calculate the mean values of the neutron energy and the yields at different angles a MATLAB function (ACEPD.M) is developed which reads from ACEFLX output. Thanks to the newly introduced target thickness condition described above, the target thickness is progressively reduced to simulate a very thin target (0.1 μ). The neutron mean energies and yields at different angles converge to the theoretical values (indicated by crosses), as shown in Figures 7 and 8. This also confirms that the random variables are properly sampled from the angular distributions.
4. The OKTAVIAN Tungsten Experiment
The OKTAVIAN facility has been devoted to perform studies on fusion neutron-related subjects by university researchers and has also served for international collaborations. The facility was a contributor to the SINBAD Project. SINBAD declared objective is the validation and benchmarking of computer codes and nuclear data used for radiation transport and shielding problems and the preservation of a unique set of experiments for the present and future needs.
The integral experiment under consideration consisted of measurements of the neutron leakage spectra from a 40 cm diameter tungsten (W) sphere pile. W is a candidate material for the first wall of fusion reactors. In the computational models, the evaluated nuclear data for W are from the IAEA files compiled in April 2007 [5]. The source neutrons were produced by bombarding a solid Ti-T target at the centre of the W pile with a 250 KeV deuteron beam. The geometry of the MCNP5 input file is based on the available information about the experimental setup. The precollimator and the main collimator bodies are modelled whereas the simulations in the original analyses neglected them. The information on the target assembly is retrieved from the report [6] and allows only an approximate model of the accelerator structures. The results presented here tend to focus on the issues arising in the analyses of MCNP5 simulations of integral experiments and directly concerning the D-T neutron source. Three approaches to the source specifications in the MCNP5 input file are assessed.
(a)In this first approach, the monoenergetic angle-dependent yields of source neutrons are given by empirical expressions which hold for thick targets. The method is illustrated by Csikai et al. [7] and is commonly used for modelling neutron sources. For comparison with measured spectra, the calculated spectra usually need to be resolution-broadened.(b)In the second approach, the use is made of the measured spectrum in the forward direction of the bare source, provided by the experimentalist. The source is defined by scaling the outgoing-energy and the magnitude of the measured spectrum such that the positions of the peak and the yield match the measured average values. In this way, resolution broadening and certain geometrical features (such as room-return, collimator, etc.) are implicitly taken into account. This procedure has been used by A. Trkov for the OKTAVIAN nickel sphere benchmark experiment [8].(c)The upgraded D-T source model described in Section 2 is applied. The deuteron energy is 250 KeV, while the other parameters are left at their standard values. An equivalent source definition (in the form of the SDEF and DS cards) is prepared for users that do not wish to recompile MCNP5 with the source routines by calculating the neutron energy distributions and yields in point detectors in void at different angles. These two forms of input give the same results and are not treated separately. The source spectra at forward angle for cases (b) and (c) are compared in Figure 9. For clarity, the plotted spectra are normalised at the peak. The scattering effects of source neutrons with the structures of the experimental setup are intrinsic to the measured response function, as can be inferred from the tail below the main peak (labelled “OKTAVIAN spectrum”). The energy spread due to the deuteron interactions inside the target as calculated by the source routine is small (labelled “SRCDX spectrum”). The calculated spectrum partly justifies the approximation in case (a), in which the energy distribution is a delta function. At other angles the source neutron distributions are similar, except for a shift in energy and the total yield.
Due to the finite resolution of the detectors, their finite size, resolution of the electronics, and so forth, some smearing in the measured spectrum is inevitable. To compare the calculated and measured spectra, the calculated ones usually have to be resolution-broadened. Assuming the applicability of the central limit theorem, the calculated spectra for cases (a) and (c) are broadened with a Gaussian resolution function. Note that resolution-broadening is not needed in case (b), because it is approximately accounted for in the source definition. The neutron leakage spectra obtained with the sources above described are compared in Figure 10. The best agreement with experimental spectrum is obtained with the recommended OKTAVIAN source. The calculated spectra at the peak agree well with the measurement. The observed discrepancy in cases (a) and (c) between 5 and 12 MeV is an indication of inadequate modelling of the detailed geometry of the experimental setup: the use of the measured response function for source definition in case (b) effectively corrects for some modelling deficiencies. Although the source in case (c) is physically better, the overall result compares worse with the measurements because it is not possible to improve the physical models of the experimental setup due to lack of information. The differences below 0.5 MeV could be due to the W nuclear data or the modelling assumptions about the surrounding structures.
5. Conclusions
The interest in the materials for nuclear fusion reactors promoted extensive research in facilities producing D-T source neutrons on small scale, as the low-voltage deuteron accelerators with Ti/T targets.
In the present work, the basic principles of the D-T fusion reaction are clarified. The calculations of the neutron energy and angular distribution are performed in the framework of relativistic kinematics. The source routines allow realistic modelling of the Ti-T target assembly, which is the starting point for the analysis of the integral experiments. The uncertainties contained in the source spectrum propagate to the neutron leakage spectrum. The analysis of the OKTAVIAN experiment showed the energy regions and the magnitude of the corresponding discrepancies which arise from the source definition and not from the nuclear data of the primary measured structural material.
In this early stage of the analysis, progress is made in establishing codes and procedures for the modelling of the D-T source neutrons. The importance of detailed and accurate specifications of the experimental setup is emphasized if integral experiments are to be used to discriminate between different sets of nuclear data.
Acknowledgments
The authors acknowledge Dr. Kodeli for providing the original source code and Dr. Pillon for his valuable suggestions. Special thanks to B. Zefran for helping with the computer environments.