Abstract
The capabilities of the SCALE6.1/MAVRIC hybrid shielding methodology (CADIS and FW-CADIS) were demonstrated when applied to a realistic deep penetration Monte Carlo (MC) shielding problem of a full-scale PWR containment model. Automatic preparation of variance reduction (VR) parameters is based on deterministic transport theory ( method) providing the space-energy importance function. The aim of this paper was to determine the neutron-gamma dose rate distributions over large portions of PWR containment with uniformly small MC uncertainties. The sources of ionizing radiation included fission neutrons and photons from the reactor and photons from the activated primary coolant. We investigated benefits and differences of FW-CADIS over CADIS methodology for the objective of the uniform MC particle density in the desired tally regions. Memory intense deterministic module was used with broad group library “v7_27n19g” opposed to the fine group library “v7_200n47g” used for final MC simulation. Compared with CADIS and with the analog MC, FW-CADIS drastically improved MC dose rate distributions. Modern shielding problems with large spatial domains require not only extensive computational resources but also understanding of the underlying physics and numerical interdependence between -MC modules. The results of the dose rates throughout the containment are presented and discussed for different volumetric adjoint sources.
1. Introduction
The Monte Carlo (MC) simulation of deep penetration shielding problems is a very challenging task. Shortcomings of the classical MC variance reduction (VR) techniques, which are absolutely necessary to get any answer at all, applied to such complex shielding problems are a much known issue. A hybrid deterministic-stochastic shielding methodology is frequently used today when facing such problems [1]. This relatively novel approach in automatic VR preparation has an ultimate goal of achieving acceptable precision of the MC results in a reasonable time. Deterministic transport theory methods [2], typically discrete ordinates , are used for numerical calculation of particle transport over space-energy domain of the problem. The solution is a mesh-based adjoint function in phase-space and it is used for source and transport biasing. Such formalism is known as CADIS [3] and it is based on the concept of the adjoint function [4] (i.e., solution of the adjoint Boltzmann equation) being the importance function toward the specified objective of user’s interest, with zero-variance solution in the limit. The typical use of CADIS is when optimization of localized results is needed, such as point or region detectors, comprising small entities (units) of much larger global unit. Generalization of this adjoint methodology for obtaining global MC distributions is far more difficult task. For the MC simulation resulting in uniformly small statistical uncertainties over large portions of phase-space overlaid with Cartesian mesh grid, additional forward calculation is needed for proper adjoint source weighting. In this case, the adjoint source is discretized over mesh cells, where every voxel has a biased strength (in space and energy) to result in global MC distribution with fairly uniform uncertainties. FW-CADIS (Forward-Weighted Consistent Adjoint Driven Importance Sampling) method [5, 6] in MAVRIC (Monaco with Automated Variance Reduction using Importance Calculations) shielding sequence of SCALE6.1 code package [7] was used with aforementioned hybrid methodology. The objective of this paper was to determine the fairly uniform dose rate distributions throughout typical PWR facility including containment building. This is representative deep penetration shielding problem, addressing modern engineering problems for mapping doses everywhere inside PWR [8–10]. Regarding the size and the model’s complexity, the utilization of the manually fine-tuned VR parameters is an impossible task. Therefore, we address the problem with the described modern hybrid shielding methodology.
The SCALE6.1 general geometry package (SGGP) was used for modeling of the typical PWR facility based on the H. B. Robinson-2 Pressure Vessel Benchmark (HBR-2) [11] critical core. Validation of the model (via detector reaction rates) was conducted in previous paper [12], where good agreement with referenced TORT results was obtained in accordance with US NRC regulations [13]. With previous HBR-2 benchmark calculations and later geometry expansion to primary loop elements [14], we now make the final geometry generalization to a full-sized, typical PWR facility with containment structure. Typical industrial and text-book data were used for dimensions and materials required: reactor internals, upper and lower reactor pressure vessel (RPV) head, biological shield, steam generators, primary pumps and pipes, concrete structures such as floors and walls, and finally containment building.
The preliminary shielding results of the PWR containment were obtained using the SCALE6.0 code but computational meshes were rather coarse and MC statistics was not satisfactory [15]. The shielding calculations in this paper are focused on hybrid capabilities in SCALE6.1/MAVRIC to produce well-converged particle fluxes and dose rates throughout containment structure originating from critical reactor core and activated coolant, which becomes additional gamma source in operating reactor. The adjoint source was examined as external air of the model and containment air to enable uniform particle attraction in all phase-space by using hybrid methodology, but the selection between CADIS and FW-CADIS resulted in drastic differences in final MC dose rates. This was notably evident for neutrons and photons originating from reactor core (compact volume) and was less evident for photons coming from activated coolant (distributed volume). The advantage of using FW-CADIS methodology with adjoint source weighting over CADIS was clearly quantified.
The paper is organized as follows. Section 2 gives the description of the SCALE6.1 code package with accent on the MAVRIC shielding sequence. Section 3 describes the PWR facility model preparation: ionizing sources, geometry, and calculational parameters. Section 4 gives MAVRIC dose rates inside PWR containment for analog MC, CADIS, and FW-CADIS methods. Mitigation of observed ray effects is considered in Section 5. Section 6 gives discussion and conclusions, while the referenced literature is given at the end of the paper.
2. The SCALE6.1 Code Package
The SCALE6.1 code package is a comprehensive modeling and simulation suite for nuclear safety analysis and design. It was developed for the US NRC for the purpose of the evaluation of nuclear facilities and radioactive package designs. The SCALE6.1 modular code is devised in analytical (control) sequences with their functional modules for performing criticality, shielding, radiation source term, spent fuel depletion/decay, reactor physics, and sensitivity analyses. For the purpose of paper clarification, only used sequences will be mentioned with the focus on the MAVRIC hybrid shielding sequence. Additional mathematical details are available elsewhere [7].
The criticality sequence (CSAS6) uses a 3D multigroup MC transport code KENO-VI to provide problem-dependent, cross section processing followed by calculation of the neutron multiplication factor . The MAVRIC hybrid shielding sequence is based on CADIS methodology and its generalization with forward flux weighting is FW-CADIS. Excellent recent papers on methodologies were done by scientists at Oak Ridge National Laboratory [16, 17]. CADIS and FW-CADIS are based on the concept of the importance function which is the solution of the adjoint Boltzmann transport equation [2, 4]. These hybrid shielding methods are used for the calculation of space-energy-dependent VR parameters in the form of the weight windows (i.e., importance map) and biased source, which work in tandem. The VR parameters are automatically transferred to functional module Monaco which is multigroup fixed-source 3D MC transport code. The integrated transport code Denovo [18] is used for the calculation of VR parameters over orthogonal meshes using Koch-Baker-Alcouffe parallel sweep algorithm and nonstationary Krylov methods to solve within-group equations. The results from even intermediate quality Denovo calculation will provide superior VR parameters compared to user’s ability to manually tune the same ones for Monaco calculations. For shielding calculations involving localized results such as point detectors or small region responses, CADIS methodology is preferred, and it is based on single Denovo adjoint transport solution. For shielding problems with multiple tallies (point and/or region detectors) or mesh tally over large portions of phase-space, an extension of CADIS method called FW-CADIS can be used to obtain uniformly small relative uncertainties [19, 20]. FW-CADIS requires additional Denovo forward transport solution to approximate multigroup fluxes and responses used for inverse weighting of the adjoint source. In case of optimizing global MC results over large meshes, every adjoint source cell becomes a discrete source element with weighted strength to ensure the same particle population for distant and close portions of phase-space. In a hybrid shielding methodology, we are interested to find a final MC solution of the steady-state transport equation:where is the linear transport operator, is the forward flux, and is the total source. A solution of (1) is based on deterministic approximation of the steady-state, multigroup adjoint transport equation:where is the linear adjoint transport operator, is the adjoint flux, and is the adjoint source. The adequate numerical solution (i.e., approximation) will provide the means to accelerate the final MC simulation via VR parameters. For global MC answers, additional forward deterministic approximation is needed to perform forward-weighting of the adjoint source. In FW-CADIS hybrid shielding methodology, the adjoint source is typically described as a product of the geometric function and the energy spectrum corresponding to the user’s function of interest (i.e., response function), which is often cross section for reaction rate calculations or dose function [7]. In that case, adjoint source weighting is done with the integral of product of the response function and the Denovo forward flux :A particle target (average) weight is used in Monaco MC simulation for particle splitting/roulette in the form of the space-energy-dependent mesh importance map and it is inversely related to Denovo adjoint flux bywhere is the normalization constant. This ensures consistency between source biasing and particle biasing, where birth weight of the particles matches target weight of importance map. Since the objective of our PWR containment calculation was to determine global dose rates inside containment structures with uniformly small uncertainties, the aforementioned FW-CADIS methodology was a highly desirable choice. To achieve such objective, one has to construct such importance function which will represent the importance of achieving uniform MC particle distribution throughout the model [9]. It was shown that this corresponds to weighting the adjoint source with the inverse of forward Denovo response. Once the adjoint source is determined, the hybrid methodology is used for the calculation of source biasing parameters and weight windows for optimized Monaco simulation. There are several multigroup cross section libraries distributed within SCALE6.1 code package. For criticality eigenvalue calculations, the “v7-238” library was used, and, for shielding calculations, we used “v7-27n19g” library for Denovo and “v7-200n47g” for Monaco. Primary data for both libraries originate from the ENDF/B-VII.0 nuclear data library [21].
3. PWR Facility Model Preparation
The simplified model of a typical PWR facility with two-loop reactor, primary loop elements, and containment building was developed using SCALE6.1/MAVRIC sequence. The reactor is the standard Westinghouse PWR type with thermal power of 2300 (710 ). The reactor core consists of 157 fuel assemblies (15 × 15 matrix, pitch 21.504 cm) and it is radially surrounded by baffle plates, core barrel, thermal shield, RPV, and biological shield. The heterogeneous fuel elements have rather complicated design which is irrelevant for the purpose of containment shielding calculations on such large scale. Therefore, the fuel elements are approximated as homogenized axial regions using referenced data from the HBR-2 benchmark [11]. The SCALE6.1/MAVRIC model of the HBR-2 nuclear reactor was validated against referenced benchmark data and was later extended to include primary loop elements [12]. The final geometry generalization now includes complete containment structure with simplified interior. Typical industrial and text-book data were used for PWR components such as reactor internals, upper and lower reactor pressure vessel (RPV) head, biological shield, steam generator, primary pumps and pipes, and massive concrete structures such as floors and walls.
3.1. Definition of Ionizing Sources
The critical reactor core was uniformly sampled in space (“flat” spatial profile) and had Watt spectrum distribution for thermal fission of 235U with a = 1.028 MeV and b = 2.249/MeV ( is the normalization constant) [7]. The total neutron intensity was n/s which corresponds to one half of total thermal power (2300 ) to better simulate neutron spatial gradient from the center to the core periphery (i.e., core self-shielding). The fission photons were also included with mean value of 7.04 per 235U fission and the secondary gamma emission in neutron transport was explicitly included in MC simulation. To include the fission photons, a mesh-based version of the neutron fission source was required and obtained using KENO-VI eigenvalue calculations. This option in SCALE6.1 is known as the CAAS (Criticality Accident Alarm System) [7] for saving space-energy fission distribution over user-defined mesh for active neutron cycles. To correctly account for how many source photons are released per source neutron, the system -parameter calculated by KENO-VI was used (2.46 neutrons/fission) with = . The activated primary coolant becomes additional gamma emitter in operating reactor, since the product of reaction with 16O in water is beta-active having a half-life of 7.13 s and additionally emitting gamma ray of high energy, 6.123 MeV, with the emission probability of 92%. The energies of the emitted photons from the coolant of a typical PWR were taken from ANSI/ANS-18.1-1999 [22], where activities for more than 50 isotopes are listed, together with plant specific scaling parameters. The details of the source term derivation can be found in the previous paper [14], together with discrete energy spectrum of the activated nitrogen which accounts for 94–97% of the total coolant activity. Using the total coolant volume of cm3, its total mass of g, and its working density of 0.79 g/cm3, we obtained total gamma source strength of photons/s. Together with photon emission spectra of , prepared for the multigroup library “v7_27n19g,” we obtained complete information for distributed coolant gamma source.
3.2. Geometry Extension
The hole option inside MAVRIC was used extensively for placing smaller units inside larger units utilizing SGGP combinatorial geometry. The homogenization process was conducted in the same way as with reactor internals, which is mass conservation based on best available input data. The MAVRIC model of the full-sized PWR reactor (coolant not shown) is shown in Figure 1 and the primary loop components are depicted without air for the purpose of clarity. Symbolically added arrows show direction of primary coolant motion in reality. Figure 1 is also showing MAVRIC model of PWR containment with front quarter removed. One can notice the reactor, the cylindrical biological shield, the primary loop elements, and the concrete shielding structures. The walls and the floors have thickness of 100 cm and 50 cm, respectively. All concrete structures were made of 02-B type concrete with density of 2.275 g/cm3. The reactor midplane is in = 0 cm plane and the primary pipes are in = 342 cm plane. The radius of containment inner steel protective shell (3.8 cm thickness) is 1600 cm and the radius of outer concrete layer was 1826 cm (76 cm thickness). The global unit, corresponding to cylindrical containment building, has diameter of about 36 m and height of about 78 m, which results in formidable MC shielding model. Although many small auxiliary components inside real containment were discarded in this phase of research, they are not relevant for the paper objective, since their inclusion would dramatically slow down MAVRIC calculations. Only gross containment structures are important because of their immense attenuating power, such as the concrete floors and the walls.
3.3. Calculational Parameters
The SCALE6.1/MAVRIC sequence was used for PWR containment shielding calculations implementing analog MC, CADIS, and FW-CADIS methodologies. We used workstation with 32 GB of RAM and Core i5 CPU. The lower limit of computer memory consumed by the Denovo state in double precision [23] can be expressed aswhere is the number of mesh cells, is the number of energy groups, is the number of Krylov vectors (10 by default), is the Legendre order of scattering cross section expansion, and is the number of unknowns per cell corresponding to spatial discretization scheme (default is Step Characteristic (SC)). The SC spatial scheme is fast with low memory requirements and robust regarding the mesh size. It always produces positive fluxes for given positive sources. However, it is 1st-order accurate compared to other high-level spatial schemes, meaning that spatial discretization error decreases linearly with mesh refinement. One can notice that Denovo model with several millions of cells and parameters with “v7_27n19g” shielding library produces state size over 30 GB of RAM, which is very memory-demanding task. Since memory consumption grows more quickly with increasing Legendre order, we gave preference to the number of Denovo mesh cells (better model voxelization) and used for all presented calculations. Accuracy of the Denovo solution in space and energy is not paramount, since only the idea about flux shape is what will accelerate final Monaco calculation. However, the more accurate the importance map is, the fewer number of histories inside MC calculation is needed to converge to the same precision, but memory consumption then becomes large. We used Denovo mesh with cells (145 × 145 × 260) and parameters and multigroup flux tolerance was set to = 10−6; the equation set was as SC, and the solver engine was GMRES method [7]. We also used the diagonal transport correction for treatment (fix-up) of zero and negative fluxes which are artifacts of finite cross sections Legendre order expansion in Denovo. The boundary conditions were all vacuum type. The implicit capture process in Monte Carlo, also known as survival biasing, was used with particle lower weight set to 10−4 for analog MC game. CADIS and the FW-CADIS were investigated with the adjoint source defined as the air inside and outside of the containment. The spectrum of the adjoint source corresponded to the user’s function of interest, which was the dose rate function: the neutron dose ID = 9029, the gamma dose ID = 9504, or the total dose ID = 9729 expressed in units (rem/h)/(particle/cm2/s). Monaco mesh had cells (112 × 112 × 200) and it was the same for all calculations. The total number of Monaco histories was (4000 batches and 25000 neutrons per bath) for all calculations. Denovo was used with “v7_27n19g” shielding library, while Monaco was used with “v7_200n47g” library to conservatively simulate low-energy particle transport (secondary gamma emission). The cross section processing codes were the BONAMI (Bondarenko factors) for unresolved resonance region and the WORKER/CENTRM/PMC (pointwise data collapsing using flux) for resonance region.
4. MAVRIC Dose Rates inside the PWR Containment
4.1. Results of the Analog MC Simulation
The total CPU time for the analog MAVRIC simulation was 1.30 days for reactor core calculations with lower particle weight set to for all 247 groups in “v7_200n47g” library. The neutron-gamma dose rates from the reactor core with relative errors () are depicted in Figure 2 in = 0 cm plane. There is practically no particle transport outside of the reactor core. These results dramatically show the necessity of VR methods in difficult shielding problems. Improvements to the analog MC simulation cannot be obtained regardless of the number of particle histories or available computer power. The same analog MC simulation was repeated for coolant gamma source, where reactor core was treated as “black” absorbing material, and the results were better in quality and quantity due to the distributed source in space (Figure 3, = 0 cm plane). The total CPU time was 1.64 days. The quality of the coolant gamma dose distribution deteriorates rapidly as the distance from the coolant increases. There is no gamma transport through the first concrete layers, so meaningful results are only in the air regions around the coolant. Obviously, the usage of the VR techniques is inevitable, especially if one seeks global MC distributions with uniformly small uncertainties.
4.2. Results of the CADIS VR Method: Internal Adjoint Source
The CADIS hybrid methodology was used to optimize MC dose rates inside air-filled regions of PWR containment, since those regions represent the most likely path of radiation streaming which is important for structural material activation and maintenance crew. The adjoint source with focus on air inside containment has a constant, uniform strength and this property is directly determined by CADIS methodology. Figure 4 is showing adjoint source in = 0 cm and = 342 cm planes for the neutron case. The spectrum of adjoint source corresponded to dose rate function of neutrons (first case) and photons (second case), respectively. The non-air-filled regions are excluded from the adjoint source and thus are presented with white color. The uniform adjoint source over air-filled regions of the containment produces fairly uniform adjoint fluxes. The inverse relationship between the adjoint flux and the particle target weight results in lack of the strong particle weights gradient in importance map. This in turn gives poor particle splitting in the regions of user interest as the distance from the forward source increases. The neutron importance map, that is, the target (average) weights for the first neutron group (6.37–20 MeV), is shown in Figure 5. Very similar distribution of target weights was obtained for the gamma case.
Since every cell of the adjoint source has the same strength, the MC particle transport exhibits self-shadowing effect. Air regions in the vicinity of the true (forward) source are getting most of the particles at the expense of distant air regions. In turn, all distant air-filled regions are shadowed by near air-filled regions exposed to higher dose rates. The neutron and gamma dose rates from the reactor core are depicted in Figures 6 and 7 with relative errors in the = 0 cm plane. The CPU time for adjoint Denovo was 2 h and it was 10 h for Monaco for the neutron case. For the gamma case, the adjoint Denovo CPU time was 2.7 h, while Monaco CPU time was about 1 day.
The independent CADIS calculation of gamma dose rates from the primary coolant was performed with the next CPU times: 26.81 min for adjoint Denovo and about 10 h for Monaco. The gamma dose rates from the coolant are depicted in Figures 8 and 9 for = 0 cm and = 342 cm planes, respectively. The better containment coverage with photons is evident for coolant source, since it is distributed over the large space.
One can notice that the importance map fails to transport particles throughout the first concrete layer in all of the CADIS cases. The Monaco results are statistically meaningful only in the surroundings of forward sources: the reactor core and the primary coolant loops. Even though CADIS is a powerful hybrid VR method, it is not created for obtaining global MC distributions over meshes but for localized answers such as point and region detectors. The reason for that is isotropic adjoint source strength in phase-space of the problem, without ability to uniformly transport particles to distant regions with many orders of attenuation. Weighting the adjoint source strength with the forward Denovo dose rates will generate space-energy-dependent volumetric adjoint source, which will be able to direct particles toward regions where their population is small. This approach requires additional Denovo calculation in forward mode for the approximation of multigroup fluxes. This leads to hybrid methodology generalization known as FW-CADIS, which is explored in the next paragraph.
4.3. Results of the FW-CADIS Method: External Adjoint Source
The adjoint source with uniform strength surrounding the containment building at the model boundary is a very attractive way for calculating dose rates everywhere inside the model. The optimized doses are expected at adjoint source location and the reasonable doses are expected in the space between. Such definition would produce importance map able to drive source particles outward from the source toward the model periphery. However, the self-shadowing effect of CADIS would offset such an approach. It would be a very difficult and iterative task for the user to manually find the right amounts of adjoint source weighting as a function of space and energy. That task becomes even harder with increasing number of cells and practically impossible if one looks for optimized values over fine-mesh with millions of cells. In order to obtain reasonable uncertainties in each cell of mesh tally covering the large model, the adjoint source in each cell should be inversely scaled with Denovo forward dose rate in that specific cell. It results in better particle transport to the low dose areas. So, with FW-CADIS, an additional forward discrete ordinates calculation is needed to approximate dose rates used for adjoint source redistribution (forward-weighting). In this way, the adjoint source is reinforced where the forward flux is low and weakened where the forward flux is high. The adjoint flux resulting from such adjoint Denovo calculation has a meaning of importance function and it is used for VR preparation. It is expected that the external adjoint source with FW-CADIS will produce better results, so in turn the CADIS external adjoint source was not investigated. Figure 10 is showing the external adjoint source in = 0 cm and = 342 cm planes optimized for gamma dose rates. It is increased in the regions with the highest gamma flux attenuation and that is predicted by Denovo to be the exterior around concrete foundations. The similar distribution in space was obtained for the neutron case. The visualization of the adjoint source in 3D was done using the program VisIt [24] and it is shown in Figure 11.
VisIt clearly indicates the complexity of the adjoint source weighing in phase-space. User’s a priori intuition and experience for manually tuning such source are futile when considering large number of mesh cells and impact of the detailed geometry on the particle transport in 3D environment. Compared to previous CADIS calculations, the importance map now has much more details regarding variable photon target weights (Figure 12, 1st group, 10–20 MeV). The global minimum of photon weights is encircling the external air near concrete foundations implying intensive gamma splitting in that region. Such weight windows will evidently not result in particle self-shadowing with first air layers in the model exterior, which was present in CADIS calculation.
Figures 13 and 14 are showing the Monaco dose rates with relative errors for neutron and photon cases in = 0 cm and = 0 cm planes, respectively. Some representative CPU times for the gamma case are as follows: forward Denovo, 2 h, adjoint Denovo, 2.5 h, and Monaco, 35.5 days. Considerable more computing effort is needed in FW-CADIS to optimize dose rates outside of the thick (76 cm) containment. The massive shielding around the reactor and the concrete floors are preventing particle transport in the foundations and the thick concrete structures around the nuclear reactor, such as biological shield and concrete walls/floors, produce heavy shielding against ionising radiation, which results with some white areas (no results). It is of interest to notice that the white areas do not represent part of the adjoint source, so the optimized results are not expected in those areas, only in the external air-filled regions. Even though reasonable results were obtained through most parts of the containment, some parts of the user’s location of interest did not receive enough particles for the valid MC statistics. This could be improved with a more thick external air layer, that is, by increasing the adjoint source volume to reinforce particle attraction in these directions.
The characteristic anomalies of the discrete ordinates algorithm used with low / parameters known as ray effects [2] are evident, especially in relative errors of the neutron case. These rays of high MC errors are rare heavy-weighted particles in Monaco coming from deterministically obtained VR parameters by low-order quadrature formula. Such formula is unable to adequately represent scalar flux from the angular one, which is calculated exactly in the discrete ordinate directions. Even with high-order solution, the flux will exhibit unphysical oscillations about a true solution, but resulting integral quantities will converge to the expected values. The containment problem is a good example, where reactor core particles stream freely through the air above the reactor head without having been scattered. Such very peaked angular distribution of fluxes is common in problems with small sources compared to large transport media that is not highly scattering. One has to note that Denovo is primarily used for quick approximation of the fluxes, not for rigorous deterministic solution, which demands cell size comparable to particle’s mean free path. In such sense, larger cells can be used and / parameters can be lower and it will still dramatically accelerate Monaco. Additional Denovo control parameters can be used to refine transport solution without burdening extra memory [14].
The independent FW-CADIS calculation of gamma dose rates from the primary coolant was performed with next CPU times: forward Denovo, 30 min, adjoint Denovo, 27.16 min, and Monaco, 17 h. The gamma transport of only first 5-6 high-energy groups, corresponding to coolant source spectrum, significantly lowers the total computational time. The gamma dose rates from the coolant are depicted in Figures 15 and 16 with relative errors for = 0 cm and = 342 cm planes, respectively. One can notice significant improvements in gamma transport through concrete structures of primary loop elements with average MC relative errors below 10%. Even though satisfactory neutron-gamma dose rates were obtained for most parts of the containment interior using weighted adjoint source at the model periphery, the weighted internal adjoint source might give more global MC results. Using the internal containment air as forward-weighted adjoint source with reciprocal of the Denovo forward dose might be a better approach, since the volume of the adjoint source is much larger. Enhancing the importance of regions with small doses while lowering ones with high doses throughout containment interior should result in the global Monaco distributions with uniformly small uncertainties. Such FW-CADIS approach is investigated in the next paragraph.
4.4. Results of the FW-CADIS Method: Internal Adjoint Source
The obtained air-focused adjoint source inside containment with FW-CADIS methodology has nonuniform strength and it is depicted in Figure 17 in = 0 cm and = 342 cm planes, respectively. Every mesh voxel represents discrete adjoint source element inversely scaled with Denovo forward total dose. The spectrum of the adjoint source corresponded to the total dose function for which the problem was optimized. This calculation implemented Monaco multisource option, where the final gamma dose is a total one, coming from the reactor and the coolant. The favorable effect of such definition could be seen in uniform particle densities inside containment, since the adjoint source strength is increasing toward model periphery in a very complicated way. For the complex MC models such as this one, it is out of the user’s ability to make any judgment on the adjoint source redistribution in phase-space. Additional understanding of underlying process can be gained with advanced computer visualization programs, such as VisIt. The 3D adjoint source depicted with VisIt in Figures 18 and 19 has a cut-plane in = 0 cm and = 0 cm, respectively, showing regular (pseudoplot) and isosurface plots. The isosurface or contour plots represent the interpolated surfaces over Denovo mesh cells with the same adjoint source strength. The portion of space with higher importance will have higher adjoint fluxes and inversely lower target weights. The complexity of the adjoint source redistribution, that is, the forward-weighting, is clearly evident. For this multisource calculation, we used slightly denser meshes for Denovo ( cells) and Monaco ( cells) with total of histories.
(a)
(b)
(a)
(b)
How the automatic procedure inside MAVRIC for obtaining deterministic solutions enables practical use of the hybrid methodology is evident and it significantly improves efficiency of shielding analyses. Such manual fine-tuning procedure would be impossible task for the user regarding the size of the model and its detailed internal structure. The visualization of such distribution reveals immense complexity which can successfully be tractable only via deterministic solution in phase-space.
The total CPU time for this multisource shielding problem was about 28 days, where Denovo consumed negligibly small part: about 3 h for forward run and 4 h for adjoint run. The bulk of the processing time was thus taken by the Monaco MC module.
The neutron dose rates from the reactor and total gamma dose rates from the reactor and the coolant with relative errors are depicted in Figures 20 and 21, respectively. The results look very promising, since most of the air-filled cells inside containment have relative errors below 20% and are comparable to the previous FW-CADIS dose rates with external adjoint source. It is advisable to increase the number of histories.
Starting with the SCALE6.1 version, several auxiliary programs became available for processing selected output files produced by MAVRIC. We used “MTSPLIT” for splitting of part of a mesh tally file (.3dmap file format) into a separate mesh tally file and then “MTBINOP” to perform binary operations on such files. We also used “MT2VTK” to convert one dataset of one family in mesh tally file to VTK formatted ASCII text file for VisIt visualization. For example, it was of interest to find the ratio of neutron and gamma dose rates originating from the reactor core. Such ratios in = 342 cm plane with relative errors are depicted in Figures 22 and 23, with a minimum in legend set to 1. One can notice the higher gamma dose rates contribution from the activated coolant in the area of the steam generator and the primary pump location, while the rest of the containment is slightly overridden with the neutron dose from the reactor.
The independent FW-CADIS calculation of gamma dose rates from the distributed primary coolant was repeated separately to investigate the effect of adjoint source redistribution. The total CPU time for the forward-adjoint Denovo was 40 min (20 min for each) and for Monaco it was 12 h. The forward-weighted adjoint source for coolant gamma dose rates using VisIt is depicted in Figures 24 and 25, respectively, with additional contour plots shown in (b). The gamma dose rates for planes = 0 cm and = 342 cm with relative errors are depicted in Figures 26 and 27, respectively. In contrast to FW-CADIS with the external air as adjoint source, the improvements in MC statistics over all of the containment interior are evident, which was the desired objective.
(a)
(b)
(a)
(b)
The FW-CADIS dose rates originating from the reactor core and from the activated coolant are depicted in Figure 28 as a function of radial distance (reactor midplane). One has to remember that the final Monaco MC results are optimized only for air regions inside containment, so some points outside the biological shield have large one-sigma bars corresponding to floor or similar concrete structures. It is evident that coolant activation mechanism plays an important role in reactor shielding design, since coolant gamma dose rates in the vicinity of steam generators and primary pumps are evidently higher (about one order of magnitude) than reactor neutron-gamma dose rates. Such behavior is clearly seen in Figures 22 and 23.
5. Hybrid Shielding Methodology and Ray Effects
Using FW-CADIS hybrid shielding methodology for geometrically small sources compared to low scattering environment will result in occurrence of rays with high MC errors, with some of them reaching values over 80%. These persistent and well-known anomalies of discrete ordinates method for particle transport are ray effects [2]. One encounters them in models with spatially small sources relative to large transport media that is not highly scattering. The ray effects in discrete ordinates Denovo are attributed to low / parameters which influence the Monaco MC simulation via VR parameters (biased source distribution and importance map). One has to remember that Denovo is not used for final transport solution but for deterministic flux approximation which will accelerate Monaco through VR parameters. It is in practice to use Denovo with lower / parameters and larger cells than in stand-alone shielding transport programs. The verification of obtained MAVRIC results was done with separate / calculation for VR parameters, but the computational Denovo mesh had to be dropped to 2.7 million cells (115 × 115 × 210), which is a factor of 2 less than a mesh with /. The Legendre order is crucial for memory consumption because of the quadratic term . In our judgment, the Denovo fine-mesh structure is more important for quality MC geometry voxelization (with macromaterial option) at the expense of lower / parameters. With macromaterial option, Denovo represents the material in each mesh voxel as a volume-weighted mixture of partial materials. Such pseudomaterial gives better representation of a cell in transport calculations involving originally different materials, such as air and concrete. Also, the ability to use Denovo with “v7_27n19g” library and Monaco with “v7_200n47g” library is an excellent way to capture detailed (i.e., conservative) contribution from low-energy particles to the total dose. The additional Denovo parameterization can benefit the final Monaco simulation on a computer with restricted memory. Some of pertinent control parameters are GMRES multigroup solver tolerance, transport correction factors for negative and zero flux fix-up, manual mesh tuning for important regions, and first-collision source option [7, 14].
The presence of ray effects inside Monaco, mapped from the Denovo calculations, was additionally investigated for the last multisource FW-CADIS calculation, optimized for total dose rates inside the containment. The Denovo mesh had about 6 million cells (with /) and the Monaco mesh had about 3.5 million cells. The normalized uncertainties (one-sigma values) from total gamma dose rates (reactor + coolant) are depicted in Figures 29 and 30 for different sigma thresholds using VisIt. One can notice the occurrence of the rays with high Monaco relative error characteristic for small quadrature order.
(a)
(b)
(a)
(b)
6. Discussion and Conclusions
The hybrid shielding capabilities and limitations of CADIS and FW-CADIS methodologies inside SCALE6.1/MAVRIC sequence were investigated on a workstation with 32 GB RAM and Core i5 CPU. We have applied the aforementioned hybrid stochastic-deterministic methodology on a realistic deep penetration shielding problem in a form of the PWR containment dose rates calculation. The sources of ionizing radiation included fission neutrons and photons from the reactor and photons from the activated primary coolant. The objective of the overall analysis was optimization of VR parameters for a fixed-source MC simulation to most effectively obtain global dose distribution over PWR containment, with nearly uniform statistical uncertainty. For that reason, the precise importance map had to be specifically tailored, resulting in a more successful Monaco simulation.
The analog MC simulation was proven to be futile and demonstrated necessity for advanced VR parameters. First, we investigated the CADIS adjoint source defined as containment air which proved unsuitable for the analysis objective, since the methodology is inherently optimized for obtaining results in a localized volume. Using CADIS for global results via mesh tally covering large spatial domain resulted in deficient Monaco dose distributions as a consequence of particle self-shadowing by air regions close to the ionizing sources. It was demonstrated that the problem was entirely dependent on successful redistribution, that is, weighting of the volumetric adjoint source, so we investigated FW-CADIS method. Using the Mesh Tally Viewer of SCALE6.1 code, one can plot the distribution of relative MC error in a mesh tally. Thus, in Figures 31, 32, and 33, we can see that fraction of the mesh tally voxels had less than a given amount of relative MC error for analog MC, internal CADIS, and internal FW-CADIS cases, respectively. One can notice extreme improvement in the mesh tally MC statistics when using FW-CADIS, so, for example, about 40% of the voxels in Figure 31 had less than 20% relative error. At the same time, CADIS and analog MC had 0% voxels with a relative error less than about 50%. The same behavior can be seen in Figure 32 for the reactor core gammas, while the primary coolant produced somewhat smoother distributions in Figure 33 as a consequence of distributed volume source.
Two approaches were considered regarding the FW-CADIS adjoint source placement: air outside and air inside of the PWR containment which were visualized with VisIt. It was demonstrated that the best (optimized) Monaco results were obtained for a particular region in question, where the adjoint source was defined, and the consistent results were obtained for in-between regions. The external adjoint source at the containment periphery produced reasonable results everywhere, but those results are optimized only for the locations outside of the containment. The volume of such adjoint source is another important factor, so by increasing it at the model boundary (i.e., more external air) the user can achieve extra particle attraction. There was initial concern that such peripheral adjoint source would sharply bias (i.e., roulette) low-energy particles in the interior model, since the probability of transporting them through thick shields is extremely low. Such particles are actually important for accounting for the total dose, since low-energy neutrons can undergo radiative capture and emit high-energy photons. However, the final Monaco results are in accordance with the results of the internal adjoint source, so that concern was not justified. The anomalies of the deterministic discrete ordinate calculation known as ray effects were explored in FW-CADIS. The mapping of ray effects from Denovo calculations to Monaco via importance map was demonstrated for the gamma case and visualized with VisIt. Additional conservatism in the final dose rates was gained by using the fine group shielding library for Monaco, which fully takes the effect of low-energy particle transport, such as () process with neutrons. Using “v7_27n19g” library for the memory demanding Denovo and “v7_200n47g” library for Monaco proved to be an excellent approach. Future releases beyond SCALE6.1 with parallel computing ability will also mitigate lengthy CPU time characteristic for these shielding calculations.
The desired MC uncertainty will dictate the optimum trade-off between the speed and the accuracy of solution, but in general it is difficult to quantify a priori the quality of deterministic solution. Using macromaterial option and other control parameters of Denovo is highly desirable to get more accurate solutions from a coarse-mesh discrete ordinates solution which will accelerate the final Monaco MC simulation.
Competing Interests
The authors declare that they have no competing interests.
Acknowledgments
This work was supported by the Croatian Science Foundation under Grant 3522.