Figure 7.7:Radial flux distribution in the test case. Coupling energy: 2.479 MeV
7.2 Estimation of high energy neutron induced reaction rates
Since most ADSs are developed for transmutation purposes, detailed burnup calcu-lations are especially important to optimize the design for MA burning. Burnup cal-culation needs the regular recalcal-culation of reaction rates to refresh the single-group collapsed cross-sections used in the Bateman equations describing the time evolution of the isotopic composition. As mentioned in Section 3.3 the high energy part of the neutron transport above the tabulated reaction cross-section is feasible only by Monte Carlo simulation. In the tabulated energy range it is straight-forward to calculate reac-tion ratesRwith the help of the track length estimator as described in Section 3.1.3 by calculating the following integral for each relevant cross-sectionsσ:
R= Z
V
Z Emax
0
σ(~r,E)Φ(~r,E)d~rdE (7.4) However, above the tabulated range the outcome of the reaction is calculated by physics models as described in Section 3.3 and the cross-sections σ are not available for the calculation of (7.4). In order to circumvent this problem a calculation scheme was developed and implemented in MCNPX with some modifications to the code.
MCNPX provides information about the residual nuclides of an interaction only through a special residual tally. This is rather suited for the simulation of thin target experiments, where the target nuclide is known and the change in the isotopic compo-sition is negligible during the irradiation. In order to make it applicable for reaction rate calculation in the residual tally a negative contribution is recorded according to
102 Chapter 7 Methods to improve the efficiency of ADS simulations
the weight of the incoming particle to the amount of the isotope which suffers an in-teraction, while the recording of the produced isotopes remains unchanged. With this method the net change in the amount of the given isotopes can be estimated, but not the actual reaction rates leading from a certain isotope to another. This problem causes some complication and needs approximations in the equation describing the burn-up, which is described below. A possible further development of the method would include the separation of residual tallies by target nuclide which would results in the separate estimation of each reaction channel.
A further problem with this kind of estimator of the isotopic transformations is that this is a fully analogue, collision based Monte Carlo estimator, which is much less efficient than the track length based estimator used in (7.4). Hence it provides results with huge relative errors especially for isotopes with low number densities since the collisions resulting in contributions to the estimation are very rare. In order to cope with this problem a special variance reduction technique, the application of so-called virtual cross sections has been introduced by the modification of MCNPX. This method is similar to the cross-section biasing decribed in Section 3.1.3 but instead of biasing a certain type of reaction, the reaction cross section of a certain nuclide is biased. This means that coefficientsciwere defined and the total macroscopic cross section used for the physics models (σit) of the ith isotope was multiplied by them: σit∗=ciσit. This correction shortens the mean free path of the particles so there will be more collisions with the ith isotope, which, of course, increases computing time. In order to preserve the reaction rates during the collision with theith isotope the particle with weightwis cut into two pieces. It takes part in the reaction and gives contribution to the residual tally with weightw1=wc
i while the other part proceeds with weightw2= cic−1
i wwithout changing its direction or its velocity. These parameters can be specified through the IDUM and RDUM cards reserved for special user options. With this method the relative errors can be reduced with several orders of magnitude for the isotopes with low number density. One has to take into account that the increase of the macroscopic cross section has the same effect as if the microscopic cross section is kept unchanged and the number density is increased (virtual number density). As the calculational time of the MCNPX is related to the number of nuclear interactions, one has to limit the virtual cross sections in order to avoid a total virtual number density orders of magnitude higher then the total number density of the material composition.
7.2.1 Application of the method for the investigation of transmuta-tion with spallatransmuta-tion source
The above described method implemented in MCNPX by the author was used by Istv´an Rovni under the supervision of Dr. S´andor Feh´er and the author to investigate the transmutation of minor actinides around an ESS-like spallation source [14]. The in-centive of the study was to investigate with contemporary calculation tools in a realistic geometry the feasibility of the direct use of spallation neutrons for transmutation since
Section 7.2 Estimation of high energy neutron induced reaction rates 103
this opportunity appeared again during the design phase of the ESS. All details of the study can be found in the conference paper [1] and [115].
As the goal was a realistic geometry, MA was considered to be placed around the ESS-like spallation source in a 30 cm thick and 60 cm high cylindrical transmutation blanket composed of fuel elements designed for 600 MWth gas cooled fast reactor (GFR)[116, 24]. This reactor type is also intended to be used for MA burning and the carbide fuel contains SiC inert matrix. The investigated MA content varied between 10%-90% while the rest of the fuel was SiC. This resulted in a kefffor the transmutation blanket of 0.17-0.96, respectively.
The reaction rates at the burnup steps was calculated by the modified MCNPX ver-sion while the evolution of the isotopopic concentration between the burnup steps was calculated by a special code developed by Istv´an Rovni. As it was mentioned above for the high energy reactions (above 20 MeV) the modified MCNPX version provided only the net production of each isotope but not the separate reaction rates. This needed the modification of the usual Bateman equation describing the change of isotopic concen-tration in the following way:
N˙i(t) =
m
∑
j=1ai jNj(t) +Fi (7.5)
or in vectorial form
N(t) =˙¯ AN¯ +F¯ (7.6)
where Ni is the number density of isotope i, m is the number of isotopes, ai j is the intensity (either decay or reaction rate) of the transformation of isotope j to isotope i and Fi is the net change (either negative or positive) per unit time of the amount of isotopeidue to high energy reactions. The general solution of (7.6) can be written in vectorial form as:
N(t) =¯ eAtN¯0+
eAt−I
A−1F¯ (7.7)
This solution can be calculated by the so-called matrix-exponential method based on Taylor-series expansion used by several depletion codes.
The application of the vector ¯F introduces an approximation since it assumes a con-stant production or consumption rate and cannot account for the change of the reaction rate during the burnup step due to the change in the quantity of the target nuclei. How-ever, the error from this approximation can be minimized by the short burnup steps. In this case the burnup step was chosen to be 0.1 year. The calculations were performed for one year of irradiation.
The calculation was performed both with the consideration of vector F and without it in order to quantify the effect of the high energy reactions. The contribution of ¯F to the total decrease of the amount of the isotopes is relevant despite the fact that only 4.74 % of the neutrons leaving the target have larger energy than 20 MeV (high energy neutrons). Comparing the change in the amount of238Pu and237Np in Fig. 7.8 for the
104 Chapter 7 Methods to improve the efficiency of ADS simulations
10 % MA case, one can observe that high energy neutrons increase the consumption of
237Np but do not change the transformation rate of the238Pu. Considering the fact that the237Np is not fissionable in the spectrum of a GFR but mainly238Pu is produced from it [24], one can conclude that the advantage of the application of a spallation source is that in this neutron spectrum even the237Np is fissionable. The same is valid for241Am and243Am, as they are also consumed in the spectrum of the spallation neutrons without the accumulation of other actinide isotopes.
Figure 7.8: Net change of the amount of the isotopes in the blanket after one year burn-up in case of 10% MA contents
However, from Fig. 7.9 one can see that the role of the high energy reaction dra-matically decreases as the MA content of the blanket and so its neutron multiplicity in-creases. These results clearly show that the direct transmutation by spallation neutrons is not economical even though the transformation of the MA isotopes is more efficient.
Although the spectrum of the fission neutrons is softer, their flux can be higher to such an extent that the overall result is higher MA consumption. However, in low multiplic-ity application the role of the high energy neutrons can be important and therefore the calculations scheme and the variance reduction method described above should be used.
Section 7.2 Estimation of high energy neutron induced reaction rates 105
Figure 7.9:Total MA consumption and the role of F vector at different MA content
Chapter 8 Summary
This thesis presents research results concerning the investigation of subcritical sys-tems by Monte Carlo methods. In the last decades a renewed interest toward the re-search on the neutronics of subcritical systems is experienced worldwide mainly due to the concept of the so-called Accelerator Driven Subcritical System (ADS). Another interest toward the physics of subcritical systems rises from the increased concern about nuclear proliferation since methods applicable for the determination of the multiplicity of subcritical systems can also be applied to determine the quantity and quality of fissile material for nuclear safeguards purposes.
The work included in the thesis covers two main field: the investigation of the neu-tron fluctuations in subcritical systems and the development of methods for the efficient simulation of an ADS. The stochastic fluctuation of the number of neutrons in the nu-clear chain reaction is the basis of the so-called neutron noise methods, which can be applied for the determination of the multiplication, the subcriticality level and other important parameters of a system.
In Chapter 4 a general algorithm is presented to provide variance reduction methods to the estimation of non-Boltzmann quantities, such as the ones describing the neutron fluctuations. The history splitting algorithm generates analogous subhistories from a non-analogous history. While the other existing methods for non-Boltzmann estima-tors are all optimized for γ-photon pulse height estimators, the main advantage of the method presented here is the applicability for higher multiplicity neutron problems and for the estimation of correlations between events.
Calculations were performed to prove the feasibility and applicability of the meth-ods. First a pulse-height-type estimator, namely a number-of-detections estimator was investigated, which is analogous to the problem of neutron multiplicity counting. It was shown that due to the calculational effort needed for the history splitting, careful optimization is needed for the efficient application of the variance reduction methods.
It was shown that the Russian roulette is inefficient in high multiplicity problems be-cause it produces large contributions at low probability. The problem can be overcome
108 Chapter 8 Summary
by using an alternative history control method, which is the limitation of the number of variance reduction events in a history. Simulation of neutron noise measurements was also performed, for which a new methodology is suggested involving the convo-lution with the source distribution. The history splitting method makes it possible to reconstruct the higher moments of the distribution, which has been demonstrated by the simulation of a Rossi-α measurement. Again, the application of the variance reduction method needs careful optimization. Surprisingly, the increased sampling of detection (detection biasing) is not efficient. Instead, the biasing of fission, which results in better sampling of the correlated particles, produces true variance reduction.
The history splitting method, which can be integrated into the general Monte Carlo program flow, has already been implemented in MCNP4C3. The application of vari-ance reduction makes it possible to calculate transport problems involving coincidence and correlation in real systems. Such calculations can be useful for the simulation of neutron noise measurement on different source driven subcritical systems. The precise simulation of such measurements can help the development of methods for the mon-itoring of the subcriticality in a future ADS. An other potential field of application is the simulation of neutron multiplicity counting systems, which are important measur-ing equipment in nuclear safeguards. Efficient calculation of full scale problems can help improve the accuracy of these systems and their capability to reconstruct the fissile material distribution in a sample.
Chapter 5 presents the results of a comprehensive set of neutron noise measurements performed at the Delphi subcritical assembly, which can serve as a validation case and a field of application for the above described methods. The measurements investigated the effect of the different source distributions (inherent spontaneous fissions and252Cf) and the radial and vertical position of the detectors. The evaluation of the measured data has been performed by the variance-to-mean ratio (VTMR, Feynman-α), the au-tocorrelation (ACF, Rossi-α) and the cross-correlation (CCF) method. Non-linear least squares fits have been applied in order to obtain the prompt decay constant (α) from the evaluated data. Theα values obtained from a singleα-fit show strong bias depending both on the detector position and on the source distribution. This is due to the pres-ence of higher modes in the system. It has been observed that the fitted α is higher when the detector is close to the boundary of the core or to the252Cf point-source. The higherα-modes have been observed also by fitting functions describing dual α-modes (α0 andα1). The fundamental mode (α0) showed much less variance in this case, but due to the insufficient time resolution the higher mode could not be determined at suf-ficient accuracy. Based on the set of measurements the α0 in Delphi can be estimated as 1704±53 s−1, whileα1 appears to be in the range of 8-9000 s−1. A successful set of measurements also gives a good basis for further theoretical investigations, including Monte Carlo simulations of the noise measurements and the calculation of theα-modes in the Delphi subcritical assembly. Preliminary results of Monte Carlo simulations per-formed by the non-analogous methodology presented in Chapter 4 are also presented.
The results prove the efficiency gain provided by the non-analogous techniques and the
109
applicability of the history splitting method for real systems.
Among the challenges of the simulation of an ADS a special case is the modelling of the neutron fluctuations in an ADS (e.g.: neutron noise measurements for reactiv-ity determination) as this requires also the description of the higher moments of the probability distributions. The high-multiplicity spallation source can have serious ef-fect on the neutron noise measurements, therefore, in Chapter 6 the physics models implemented in the MCNPX code have been used for the investigation of the energy correlation in a spallation neutron source. After a successful validation of the model with a neutron multiplicity measurement, the one- and two-particle energy distributions and spectra of the spallation neutrons have been determined. It has been shown that the one-particle energy distributions highly depend on the number of neutrons produced in a source event. Fewer neutrons result in a higher probability of high energy neutrons as the excitation energy per neutron is higher. It has also been observed that the averag-ing effect of the large number of secondary reactions in thick targets reduces this effect.
The energy correlation between the neutrons from low multiplicity spallation events has been demonstrated also by the two-energy distributions. Finally, one can conclude that the energy correlation is an important effect in a thin spallation target and for the com-plete description of the neutron fluctuations, the two-energy spectrum χ(E1,E2)has to be used. On the other hand, this effect reduces in thick targets and the factorization of the two-energy spectrum χ(E1,E2) =χ(E1)χ(E2) appears an acceptable approxima-tion. A methodology is presented to quantify the actual effect of the observed energy correlations on the higher moments of the neutron fluctuations in a subcritical system.
The development of coupled Monte Carlo and deterministic calculation methods has outstanding importance in the research of the ADSs, because fast and flexible tools are needed for design calculations. In Chapter 7, the methodology of an energy based volumetric and angular coupling is described and its feasibility is shown with test cal-culations. Due to the much faster deterministic calculation in the low-energy region, huge computer time can be saved. The presented calculational method can be easily in-tegrated into a burnup calculation scheme. With the proper choice of the deterministic code, dynamics calculations are feasible, as well. An other possible field of applica-tion could be the shielding calculaapplica-tions for high-energy facilities (accelerator, spallaapplica-tion source) because the deterministic codes are more efficient in deep penetration problems.
Furthermore, a variance reduction method based on virtual cross-sections has also been presented to improve the efficiency of the estimation of the high energy neutron induced reaction rate. It was shown that the application of such methodology is unavoidable in cases when the high energy spallation neutron has a significant contribution to the reac-tion rates.
The new research results of the author presented in this thesis can be summarized in the following propositions:
1. I was the first to develop a method which makes it possible to use non-analogous techniques in the Monte Carlo estimation of correlations between detection events
110 Chapter 8 Summary
and therefore in the simulation of neutron noise measurements. The so-called history splitting reconstructs all the possible analogous subhistories from a non-analogous particle history in the form of a subhistory matrix. Based on this sub-history matrix the contributions to a correlation estimator can be calculated. [5, 2]
2. I was the first to demonstrate that the use of Russian roulette is problematic when the history splitting or a similar technique is applied for the non-analogous Monte Carlo estimation of a non-Boltzmann quantity, because it results in low proba-bility, high importance contributions. Alternatively, I have suggested the use of alternative history control techniques, namely the limitation of variance reduction events in a history and in a particle track. [4, 5, 2]
3. I have suggested a completely new method for the Monte Carlo simulation of the results of neutron noise measurements based on the convolution of subhis-tory contributions with the source distribution. This method largely improves the statistics of the estimation even for analogous simulations. I have implemented the method and tested with the autocorrelation technique.[2]
4. I was the first to implement the methods in 1-3. for neutron transport in the three dimensional, continuous energy, general purpose Monte Carlo code MCNP, which made it possible to perform non-analogous Monte Carlo neutron noise simulations for real systems. I have demonstrated the applicability and efficiency of the methods with simulations for the neutron noise methods performed on the Delphi subcritical assembly. [4, 6]
5. I have demonstrated the spatial dependence of the higherα-modes in a thermal system by neutron noise measurements at the Delphi subcritical assembly. This set of measurements provides a good basis for further theoretical investigations of theα-modes.[6]
6. With the help of the physics models implemented in the high-energy Monte Carlo particle transport code MCNPX, I was the first to demonstrate that the energies of neutrons produced in spallation reactions from a single proton cannot be con-sidered as independent. The energy correlation between the neutrons diminish as the target thickness and the number of secondary reactions increases. [8, 9]
7. I was the first to develop a volumetric and angular coupling between high-energy Monte Carlo and discrete ordinates transport at the energy boundary of the high energy physical models and the tabulated cross-sections. This makes it possible to apply the computationally expensive Monte Carlo transport only at the high energies for the calculation of the spallation source and perform the calculation for the source driven subcritical core with a more efficient deterministic method.
[10]
8. I was the first to develop a methodology and a variance reduction method based on virtual cross-sections for the efficient estimation of reaction rates induced by
111
high energy spallation neutrons and implement it in the high-energy Monte Carlo particle transport code MCNPX. The method was proved to be applicable for problems, where high energy particles have significant contribution to the reac-tion rates. [1]
Acknowledgements
It is my pleasure to express my gratitude to everyone for the support, encouragement, pieces of advice and all the other help I have received during my professional carrier, and without which this PhD thesis could not be completed.
I am grateful to Prof. Dr. Gyula Csom, who had been my supervisor since my undergraduate years and who directed my interest towards the topic of transmutation and so to the Accelerator Driven Systems.
I would like to express my gratitude to all my colleagues at the Institute of Nuclear Techniques for the friendly, helpful and professional working environment in which I had the opportunity to work during these years. Especially to Dr. S´andor Feh´er, who has introduced me to the world of the Monte Carlo method, and helped me with his valuable advices. My first stay in Delft would have never been realized without his efforts. To Dr. Szabolcs Czifrus, who was always available for the language review of my English writings (such as this very thesis) and who kindly accepted to be an opponent at the internal defence of this thesis. Together with Dr. G´abor P´or, the other opponent, they acted almost as my advisors and their valuable remarks improved the quality of this final version.
I owe my gratitude to Dr. Jan Leen Kloosterman, who made it possible for me to stay at the Delft University of Technology and trusted my first, vague (and, as turned out, somewhat flawed) ideas about variance reduction for neutron noise simulations.
This initial trust was essential for the realization of all the results included in Chapter 4.
His support and enthusiasm was also a prerequisite for the measurements at Delphi.
Together with him, I am also grateful to all the colleagues I met in Delft and who made me those months a pleasant memory.
I would like to thank Dr. D´avid L´egr´ady the pool and squash games and the fruitful conversations about the mysteries of non-Boltzmann estimators both in Delft and in Budapest.
I am grateful to Dr. Imre P´azsit for his kind attention and the fascinating scientific discussions. Especially the one during a coffee break of the Physor 2010 conference, which initiated the research presented in Chapter 6. I am also grateful to all the