Radiation transport out from the reactor core: to decouple or not to decouple

. In the framework of the extension of the lifetime of currently operating reactors as well as of issues connected to decommissioning, accurate calculations of neutron and gamma responses outside the reactor core are increasingly being sought. Recently Monte Carlo calculations have been extended to provide a deep penetration capability incorporated within the eigenvalue calculation. This allows, in principle, neutron and gamma ray responses quite far outside the fissile region to be calculated within the same source-iteration scheme employed to define the neutronic responses in the fissile zone. In this paper, the new algorithm is compared to the classic decoupled approach - an eigenvalue calculation followed by a fixed source one - with the point of decoupling chosen as the fission sites. Two contrasting sample problems are discussed: a small fast research reactor and a large GEN-III Pressurized Water Reactor. The latter problem highlights the role of superhistories in maintaining the fundamental mode.


Introduction
In the last years work has been underway to extend the lifetime of currently operating reactors. In the case of Pressurized Water Reactors (PWR's), the radiation damage to the pressure vessel is a prime example of such evaluations. Furthermore there are an increasing number of requests to evaluate a variety of ex-core radiation responses linked to added safety devices or appraisals (ex-core neutron detectors, concrete ageing, detectors employed in severe accident scenarios, etc.). As far as decommissioning is concerned, accurate assessments of ex-core activation are required to develop dismantling and waste storage strategies, etc.
Normally the calculational approach to such problems, in particular when employing Monte Carlo, is to decouple the problem into an eigenvalue calculation of the fissile zone and a fixed source calculation going outside the core. The point of decoupling varies: it may be the flux leakage from the fissile zone or the fission sites within the fissile zone. Of these two options, employing the fission sites allows the number of phase space variables to be reduced to three: viz. x,y,z. (In particular the delicate question of discretization of the angle is avoided as fission is isotropic.) Notwithstanding, the fission site option still has several disadvantages: -the spatial discretization still introduces some approximations that are difficult to quantify; -for a core containing a variety of fuels and enrichments (and for a burnt core), it may not be feasible to tally the fission source at a level that is consistent with the material geometry description; -to the authors' knowledge, current public Monte Carlo codes do not allow a complete pin-by-pin description of the fission source -any user-written modifications reduce the quality assurance; -the whole process (eigenvalue calculation to generate the fission sites, interface code to create the source for the fixed source calculation, the fixed source calculation) is cumbersome and long, both in terms of computer and of human time.
Recently a new algorithm has been proposed to perform variance reduction (VR) within the sourceiteration algorithm and therefore to allow in principle excore radiation responses to be estimated within a single eigenvalue calculation ( [1][2][3]). In [1] and [3] the feasibility of the new approach was demonstrated in that localized in-and ex-core responses were calculated with a low error within the eigenvalue calculation. However the sample problems contained just a single ex-core problem (the rest were in-core).
Work is currently underway to test the new algorithm on ex-core problems in a more systematic fashion. In this paper the new algorithm is compared with the classic decoupled approach employing the fission sites at the point of decoupling for two contrasting sample problems: a small fast research reactor and a large GEN-III thermal power reactor.

The Direct Statistical Approach
The Direct Statistical Approach (DSA) was originally developed for fixed source problems [4]. It can be distinguished from classical VR approaches such as the weight window generator [5] by the specific characteristic of estimating the second-moment and therefore the error.
The VR technique that the DSA aims to optimize is splitting and Russian roulette (RR) at surfaces in space and in energy (i.e. a 2-dimensional function), either independent of, or dependent on, the weight of the progenitor arriving at the surface [4]. Because the DSA uses the second-moment, in contrast to classical approaches that rely on the first-moment, it requires also the time. This is so as to be able to find the optimum trade-off between over-splitting that takes too much time (and therefore does not run a sufficient number of histories) and under-splitting that has too high a variance per source history.

Applying the DSA to the source-iteration algorithm in eigenvalue problems
The key aspect in applying the DSA to the sourceiteration algorithm is the multi-response capability [4]. In applying VR to the eigenvalue problem, we are changing the sampling in the fissile zone and in the reflector and thereby distorting the fundamental mode. But the fundamental mode is the source for the next cycle in which again we are changing the sampling to improve the statistics of local detectors, be they in-or ex-core and logging the fundamental mode, the source for the next cycle, etc. Clearly there is a trade-off between the amount of distortion of the fundamental mode (to be minimized) and the alteration of the number of tracks at the local detector (to be maximized).
Thus we must add responses in the fissile zone (normally the fission neutron production rates in a sufficient number of spatial segments to mock-up the fundamental mode) which we name the global responses, to the local in-or ex-core responses for which VR is required. The DSA aims to optimize all these responses simultaneously.
It may be mentioned that first-moment approaches to optimizing VR can also contain a multi-response capability [6,7]. In principle therefore such approaches could be employed in the same way in eigenvalue problems. Notwithstanding, comparisons between the DSA and classical approaches have only been made so far for fixed source problems [8] apart from a single assessment for an eigenvalue problem in [1].

Superhistories
Superhistories [9] lengthen the number of fission generations in each cycle, a cycle being defined as between normalization points, in other words an "iteration" in the source-iteration scheme. Note that the usual source-iteration scheme has one fission generation per normalization cycle. Although superhistories were introduced in [9] principally to mitigate the problem of underestimation of errors in the source-iteration scheme (and to a lesser extent the problem of response underestimation), in [1][2][3] they were reintroduced to improve the variances of local responses. Furthermore the VR parameters were varied according to the fission generation within the superhistory. This was successful for in-core responses, although the gain in efficiency was only a factor of 2 or so. Instead for ex-core responses there was virtually no gain. This was expected as it is the fissions in the core periphery that contribute to ex-core responses. Thus neutron transport within the fissile zone is unimportant for the local ex-core responses.
Notwithstanding a remark in the conclusion of [1] -"Allowing a number of fission generations between normalizations of the fission source (i.e. superhistories) may help in maintaining the fundamental mode ……"the working hypothesis from that work was that superhistories of a single fission generation were appropriate for ex-core local responses whilst for in-core local responses, more than one fission generation per superhistory may improve somewhat the efficiency.
In the following sample problems, a single fission generation per superhistory was first chosen. This choice was satisfactory for the TAPIRO problem with a small, fast, core. The benchmark results of this problem provide the first main result of this paper. Instead for the large GEN-III thermal core, this choice produced results which were flawed. The second main result of this paper is the analysis of the cause of such flaws and subsequent indications of how to eliminate them.

TAPIRO fast reactor sample problem
TAPIRO is a helium-cooled 5 kW fast research reactor at the ENEA Casaccia centre, near Rome. It has a small, cylindrical (approximately 11 cm high, 12.5 cm diameter) highly-enriched core. Because of its relatively hard spectrum, it has recently been proposed to employ TAPIRO to irradiate actinides and measure their capture and fission cross-sections. To this end, a characterization was made of three experimental channels: diametral, tangential and radial, from the fissile region through the copper reflector, up to the concrete shield (see Figs. 1-3) [10]. In each mesh interval in the filled channels in Figs. 1-3, the total capture and fission reaction rates for the following isotopes were required: 235 U, 238  Four calculational approaches were followed: 1) A single "brute force" analog eigenvalue calculation, using MCNP6.1 [11].
2) The classic decoupled approach employing the fission sites as the point of decoupling. For the second, fixed source, part of the calculation, analog VR was employed. Both parts of the calculation used MCNP6.     Further aspects of the calculational approaches 2) -4) are described as follows: Approach 2): The spatial binning of the fissile zone was 2-dimensional, r-z, with 19 axial and 14 radial bins. This could be exactly modelled with the SDEF source card in MCNP.
22 reaction rates (the capture and fission rates of the 11 isotopes) are requested in each of 226 meshes (96 meshes in the diametrical channel, 90 meshes in the tangential channel and 40 meshes in the radial channel). This is too large a number of responses for the DSA [4] which, in contrast to first moment approaches, must  240 Pu and 244 Cm and fission in 243 Am) and the response meshes were merged to make a total of 12 macromeshes (4 in each channel). Thus the problem was reduced to finding the optimum VR parameters in 16u10 phase space cells in the presence of 5u12 responses. As no weight variations were present other than those due to phase space surface splitting/RR, the weight-independent DSA was employed (see for example [13]). The resulting importances in space/energy were converted to a weight window and run with standard MCNP6.1.
Approach 4): In this case there is of course no source binning. The same 16 spatial regions and 10 neutron energy groups were employed as in the second part of the decoupled case in approach 3). Furthermore the same 5u12 responses were used as in the decoupled case. However now the fissile zone was divided into 8 segments and the fission rate was scored in each segment so as to maintain the fundamental mode. These 8 responses were added to the 60 ex-core responses. Again like the decoupled case, no weight variations were present. A superhistory of a single fission generation was employed.
In Figs. 5-7 are shown the results for the capture rates in 238 U, 242 Pu and 244 Cm respectively in the upper part of the diametral channel (see Fig. 1) for the four approaches: single eigenvalue analog, decoupled analog, decoupled VR and single eigenvalue VR. In Fig. 8 are shown the results for the capture rates in 240 Pu in the radial channel (see Figs. 1 and 3) again for the four approaches. Instead in Fig. 9 we see the fission rates in 245 Cm in the tangential channel (see Fig. 2).
Although not part of this comparison, it is of interest to note that the computer times were: 548.5 h, 638.9 h and 1071.6 h for the decoupled analog, decoupled VR and single eigenvalue VR calculations respectively. (The times with VR do not include the times required to generate the VR parameters which were however relatively small. The same does not apply to the human time of course.)     For the small TAPIRO core and given the relatively fine source binning, the approximation involved in the decoupling can be neglected. Thus the results look to be at an appropriate benchmark level for the single calculation approach without decoupling for this kind of problem.

Gen-III power PWR
Axial and radial sections of the simplified PWR GEN3 model adopted with a heavy steel reflector are shown in Figs. 10-13. The detectors were situated outside the pressure vessel radially, between the vessel and the concrete and at various depths within the concrete, and axially, above the concrete at the bottom of the well and embedded in the concrete basemat. As well as neutron responses, gamma dose rates were also requested.
A pin-by-pin description of the fuel assemblies was adopted. Given that the core geometry is defined pin-bypin, we wished also to define the fission source at the same level. Thus in contrast to the TAPIRO problem, the mechanism of decoupling of the calculation using the fission sites was relatively complex.
We did not attempt to run an analog calculation to calculate the fluxes and responses as the penetration was too great, especially for the responses embedded in the concrete basemat at the bottom of the well.

The decoupled calculation: writing the source
The core contains 241 fuel assembles, each assembly with 17u17 fuel pin positions. After taking away the guide tubes (24 per assembly), we have a total of 63865 fuel pins (including Gd-doped ones). The fissile zone is approximately 4.2m high. Each pin must be divided into a sufficient number of axial segments to model reasonably the axial power profile. We used 37 axial segments, of variable length, but the same for all the pins.
Thus we had a total of 63865u37=2363005 tallies for the neutron source. Because of memory issues, it was necessary to divide the first part of the decoupled calculation that writes the neutron source into six parts, each part estimating the fission neutron source in a subset of the fuel pins. Each of these six calculations were identical (apart from the tallies), and starting from a converged fundamental mode fission distribution, ran 4000 cycles, each of 400000 source fission neutrons.

The decoupled calculation: reading the source
To read the resulting fission source distribution, MCNP does not allow a definition of the source at the fuel pin level when more than one fuel assembly is present. Thus a user-written source was required. The authors actually employed a dummy standard source and modified the code. Rather than keep all source date in memory, for each source history the spatial coordinates were read from a direct access file.   Given that both neutron and gamma total fluxes, flux spectra and responses were requested at a variety of positions in the pressure vessel well and in the concrete surrounding the well, it was decided to divide the second part of the decoupled calculational suite into two: a neutron-only run and a n-gamma run with only gamma responses. In both cases the DSA [4] was employed. To generate the VR parameters, the number of neutron responses was reduced to 48 (with 44 spatial regions, 8 of which covered the fissile zone, and 7 energy groups with limits: 5, 2, 0.2, 0.02, 0.001 and 1E-6 MeV) and the number of gamma responses was reduced to 22 (with the same 44 spatial regions and 7 neutron energy groups plus 6 gamma energy groups with limits: 5, 3, 1.5, 0.7 and 0.2 MeV).
Having generated optimum or near-optimum VR parameters in the neutron-only and n-gamma cases, runs were made with these parameters, employing 390u10 6 source histories in the neutron-only case and 100u10 6 source histories in the n-gamma case. The results for all the requested responses, nearer or farther from the pressure vessel, were satisfactory and within the project constraints. At that point it was necessary to evaluate the uncertainties in the results arising from the decoupling.

The single eigenvalue calculation
It was desired to evaluate the approximations of a single radial bin for each fuel pin (which should underestimate the ex-core responses) and of the axial discretization of the source into 37 bins, by making two single eigenvalue calculations, one for the neutron-only case and the other for the n-gamma case.
The optimization of the DSA neutron and n-gamma calculations closely followed the fixed source cases and will not be described in detail. A key difference from the fixed source cases was the addition of 12 fission heating responses over the fissile zone making a total of 60 responses in the neutron-only case and 34 responses in the n-gamma case. This meant that the core spatial division into 12 segments to provide responses to maintain the fundamental mode was not the same as the core spatial division (8 segments) to control the neutron transport within the core. Although probably not critical, this is undesirable and was corrected in the later calculations. For both the neutron-only and the n-gamma cases the spatial regions and energy groups employed for VR were the same as in the decoupled cases. Again a superhistory of a single fission generation was employed.
Having reached the fundamental mode, cycle lengths of 350000 and 150000 fission neutrons were used in the neutron-only and n-gamma cases respectively.

Comparison of the results of the decoupled and single eigenvalue calculations
In Figs. 14 and 15 the axial profiles of the total neutron and total gamma fluxes radially just outside the pressure vessel are compared for the decoupled case and for two statistical samples (i.e. a different number of cycles of the same calculation) of the single eigenvalue case.
We see in these figures, especially in the neutrononly calculation, that there is a large distortion of the axial profile that is swamping the effect that we are looking for (the approximation of the decoupled calculation). In the rest of this section, we shall analyse the cause of such distortion. The analysis will focus on the single eigenvalue case.

Analysis of the single eigenvalue results
We suspected that the fundamental mode was not being properly maintained in the case of the single eigenvalue calculation. To examine this hypothesis, rather than look at the fluxes outside the vessel, we investigated the leakage neutron flux from the fissile zone into the radial reflector (see Figs. 12 and 13). We firstly studied the evolution of the axial profile of the leakage neutron flux with the number of cycles run for two cases: -a standard MCNP calculation with analog VR parameters and cycle size 350000 fission neutrons.
-a single eigenvalue calculation with a superhistory of one fission generation and a cycle size of 350000 fission neutrons. The same ex-core spatial cell configuration, energy group structure and ex-core responses were maintained from the previous calculations. Now instead the in-core spatial division and global responses were changed so that an identical spatial division of the fissile zone provided responses to maintain the fundamental mode as controlled the neutron

ICRS-13 & RPSD-2016
5007 transport within the core. This division was: three radial segments: the outer assembly ring, the next two assembly rings and the inner assemblies and five axial segments: < -147, -147 --63, -63 -+63, +63 -+147 and > +147 (all cm with respect to the core mid-plane). Thus there were 15 core regions and 15 global, core responses making a total of 63 responses. There were correspondingly 51 spatial regions and, as previously, 7 neutron energy groups. As before the VR parameters were then optimized to the neutron local ex-core responses and global in-core responses simultaneously.
A superhistory of a single fission generation was employed.
In Fig. 16 we see the evolution of the axial profile of the core neutron leakage flux for eight cases between 150 and 1500 cycles, with standard MCNP and analog VR. The axial profiles from the DSA run with exactly the same number of source fission neutrons are shown in Fig. 17. We see that unlike standard MCNP and analog VR. the DSA does not sufficiently maintain the fundamental mode.  Possible solutions to improve the convergence of the fundamental mode with the DSA were considered: 1) Increasing the number of in-core responses. Currently 15, the in-core responses might not have a sufficient weight in the minimization compared with the 48 ex-core responses.
2) Increasing the length of each superhistory: this should automatically raise the weight attached to the incore responses due to transmission through the fissions of the response importances in future fission generations within a superhistory. It may also allow a better settling of the fission distribution due to longer in-core transport between normalization points.
The second option was adopted and superhistories of 1, 4, 12 and 25 fission generations were run. For each superhistory case, VR parameters were generated with the DSA in the 51u7 phase space regions including the 15u7 phase space regions in the core. Given the k eff of the system of around 1.0132, fission generation numbers of: 150, 300, 450, 600, 750, 900, 1200 and 1500, were mocked up as normalization cycle numbers of: For superhistory lengths > 1 these normalization cycle numbers would only approximate the required numbers of fission generations for analog importances. It was actually not possible to maintain the required numbers of fission generations because of splitting/RR between space/energy regions in the fissile zone. Furthermore there were many more ex-core collisions with the DSA compared to MCNP analog because of optimizing to the ex-core detectors.
There were 350000 fission neutrons starting each cycle in all four superhistory cases.
In Fig. 18 we see the evolution between 11 and 108 cycles of the axial profile of the core neutron leakage flux with the DSA for a superhistory of 12 fission generations, with again 350000 source fission neutrons per cycle. We see in Fig. 18 that the convergence of the fundamental mode has markedly improved. We now analyse further the effect of the length of the superhistory.

The impact of the length of the superhistories on the convergence of the fundamental mode
Figs. 19 show the axial profiles of the radial leakage fluxes of MCNP analog and the four DSA superhistory cases for cycle numbers of: 150, 300, 450, 600, 750, 900, 1200 and 1500 respectively, or their equivalent as discussed above for superhistory lengths > 1. The results are given as a ratio to the results of a separate, independent, analog MCNP run of 3500 cycles of 350000 neutron each. We see again that a superhistory of one fission generation does not sufficiently maintain the fundamental mode. From these figures we may infer that whilst a superhistory length of 4 is markedly better than 1, a length of 10 or so should be sufficient.
To try to establish a more level playing field for these results, we show in Table 1 a number of standard counters for the final calculations in Fig. 19h. The "loss to fission" row gives some measure of the statistical effort to maintain the fundamental mode. Although the DSA with a superhistory length of 1 has the lowest value (one quarter of analog MCNP), it is clear from Figs. 19 that four times the effort will still not bring these results near the quality of the others (and would anyway take four times the time of the other superhistory cases).         To illustrate further the effect of superhistories in better maintaining the fundamental mode when variance reduction is present, we show in Figs. 20 the phase space cell importances (the z axis) found with the DSA for the four cases of superhistory lengths: 1, 4, 12 and 25 fission generations. The fifteen core regions, defined previously, have five axial segments ("z" in cm, going from bottom of core at NE to top of core at SW in Figs. 20) and three radial ones ("Inner" assemblies at NW, "Middle" assemblies and "Outer" assemblies at SE in Figs. 20). We remind ourselves that the seven neutron energy groups have limits: 5, 2, 0.2, 0.02, 0.001 and 1E-6 (MeV). Thus in Figs. 20 for each of the 5u3 core regions, there are 7 histograms, going from the highest energy red to the lowest energy blue.      Fig. 20a have higher importances at low energies of the external assemblies at the bottom of the core compared with the top. This reflects the location of the deeper penetration ex-core local responses at the bottom of the well in Fig. 11.

Conclusion
The introduction of variance reduction within the Monte Carlo source iteration algorithm opened the path to the possibility of calculating ex-core responses at deep penetrations without decoupling the calculation with the inherent approximations associated with the decoupling. Whilst for in-core local responses it had been shown that employing superhistories [9] and varying the VR parameters according to the fission generation within the superhistory can improve the multi-response variance, for ex-core responses no such gain had been found [1][2][3]. This was understandable because only peripheral fuel pins contribute to ex-core responses. Thus it had been imagined that superhistories need not be employed in excore problems and for in-core problems their only role is to improve the efficiency.
In this paper two sample problems involving ex-core responses are considered: the first a fast research reactor with a small highly enriched core, the second a simplified model of a large GEN-III PWR.
The first sample problem showed excellent agreement of the ex-core responses between the single eigenvalue run with VR (and a superhistory of a single fission generation), the decoupled run with VR and an analog decoupled run.
The second sample problem instead showed serious large-scale discrepancies between the single eigenvalue run with VR (and a superhistory of a single fission generation) and a decoupled run with VR. These discrepancies were not due to the approximations inherent in the decoupling but were caused by a flaw in the single eigenvalue calculation that did not maintain the fundamental mode sufficiently well. It was shown to a reasonably high probability that this was due to the use of the single fission generation per superhistory. Just four fission generations per superhistory improved greatly the stability of the fundamental mode, whilst 12 and 25 fission generations per superhistory improved it further. Probably 10 fission generations per superhistory would be sufficient (which is the same number recommended in [9] for other reasons).
Although only ex-core responses were considered here, the question of stability of the fundamental mode is of equal, if not greater, importance for in-core responses. Thus when employing variance reduction with the DSA within source-iteration eigenvalue calculations, superhistories with 10 or more fission generations should always be used for in-and for ex-core responses. With this caveat, a single calculation for ex-core responses should be more accurate and faster (in terms of human time) than a decoupled calculation, as well as involving a higher level of quality assurance.