Fast three-dimensional radiation field characterization in large-scale nuclear installations

Fast three-dimensional radiation field characterization in large-scale nuclear installations 1 corresp first-author Corresponding author, yuanjiebi@gmail.com Corresponding author, yuanjiebi@gmail.com yuanjiebi@gmail.com 1 1 1 1 1 China Institute of Atomic Energy, Beijing 102413, China China Institute of Atomic Energy, Beijing 102413, China The requirement of the fast three-dimensional radiation field calculation is raised during the decommissioning of large-scale nuclear installations. The most often used methods, such as the Monte Carlo and the discrete ordinates methods, have shortcomings in their simulations of such problems. The coupled Monte Carlo-Point Kernel computational scheme is developed to meet the requirement. The facility is separated into the source region and the bulk shielding region, with a common interface. The transport within the source region is simulated using the Monte Carlo method, which is by nature good at treating complex geometries. The radiation field in the bulk shielding region is treated by the point kernel approach to avoid the extremely expensive computation for deep penetration problems. The flow rate through the interface, which is given by the Monte Carlo simulation, is considered as the equivalent surface source for the point kernel calculation. Test calculations from simplified typical waste storage facilities have been performed to validate the coupled scheme by comparing them with the results conducted by the Monte Carlo method directly. The good agreement of the results, as well as the significant saving in computing time, indicates that the coupled method is suitable for the fast three-dimensional radiation field calculation. Fast radiation field calculation Point kernel approach Monte Carlo method Coupled computational scheme 1 Introduction 1 During the decommissioning of nuclear installations, the shielding conditions and the sources change quickly with the dismantling of the devices. The disappearance of one risk always makes another one become progressively dominant [ 1 1 ]. Meanwhile, measurements cannot always be done as desired. Therefore, fast simulation in a changing radiation field is quite necessary to follow up the work progress. Furthermore, traditional dose assessment methods, which often give the dose rate in some representative positions, cannot reveal the dose variations in the full space. The given dose rate may not indicate the actual value for a worker standing in a specific location. Given the situation, the need for full three-dimensional (3D) radiation field simulations is raised. The fundamental tool to solve radiation transport problems is the Boltzmann equation, which is, in general, extremely difficult to solve. The deterministic discrete ordinates method and the stochastic Monte Carlo method are the most often used methods to solve this equation [ 2 2 , 3 3 ]. The discrete ordinates method, which solves the Boltzmann transport equation by the discretization of the angular variable, could deal with deep penetration problems but has limitation on treating complex geometries and sources. The Monte Carlo method is naturally suitable for the complex geometry and complicated source modeling, however, it requires great computational expense for large- scale bulk shield, even in the cases where different kinds of variance reduction techniques are used. The discrete ordinates adjoint function could also be used to do automated variance reduction of Monte Carlo shielding calculation, which needs the discrete ordinate code to work first [ 4 4 , 5 5 ]. To meet the requirement for fast calculations, some tools based on the conventional point kernel method are also developed [ 1 1 , 6 6 - 10 10 6-10 ]. These tools are supportive of the ALARA (As Low As Reasonably Achievable) approach. The point kernel approach is fast enough even for 3D interactive simulations as compared to the other two methods. However, large numbers of approximations make the calculation accuracy very low, especially for complex geometry configurations and for the detectors close to sources. Given the difficulty in dealing with the large-scale bulk shield simultaneously while treating the complex geometries and sources using a single method, the coupled computation scheme has been developed. Dividing the geometry into a complex source region and a bulk shield region, the radiation field could be solved by combing two different methods. The complex source region is simulated first to get the particle information on the interface, which is then used as the source term of the bulk shield region. The coupled Monte Carlo-Discrete Ordinates computational scheme, which uses the Monte Carlo method to deal with the source region, and the discrete ordinates method to treat the bulk shield region have been proposed in some papers [ 11 11 - 14 14 11-14 ]. However, this coupled method is not applicable to do the fast simulation, since discrete ordinates calculations are usually computationally expensive and need skills to speed up the convergence of the iterative solution [ 15 15 ]. In view of the situation described above, we developed a fast 3D radiation field simulation method based on the coupled Monte Carlo-Point Kernel computational scheme. It could make the 3D simulation much more accurate than that given by the simple point kernel method, and it is faster than both Monte Carlo simulation and Monte Carlo-Discrete Ordinates calculations. The idea of the coupled computation scheme was initially validated by the simulation of a single radioactive waste storage container in the previous study [ 16 16 ]. In order to adapt it to much more general situations, the physical models and the algorithms, as well as the key computing modules, have been improved in this paper. Two representative examples, covering the symmetric and asymmetric sources and geometries, have been examined to validate the program system. 2 The coupled Monte Carlo-Point Kernel program system 1 As a new way to obtain the 3D radiation field, the coupled Monte Carlo-Point Kernel computational scheme is characterized by the fast and full 3D calculation. An interface is selected to separate the whole geometry into two parts: the Monte Carlo region with complex sources and geometry, and the point kernel region with simple geometry and bulk shield. An interface program is developed to convert the Monte Carlo output into the equivalent surface source of the point kernel calculation. Combing the Monte Carlo and point kernel results, the 3D radiation field of the facility can be obtained. There are six steps to do the coupled Monte Carlo-Point Kernel simulation: • Select the interface; • Do the Monte Carlo simulation; • Set up the interface program; • Do the point kernel approach; • Consider the wall/floor reflection; • Combine Monte Carlo and point kernel results to get the 3D radiation field. 2.1 Selecting the interface 2 Since the symmetry is usually present in practical problems, the solid boundary is often selected as the interface, which could significantly simplify the calculation. For the case of spherical, cylindrical, rectangular, and conic containers filled with radioactive nuclides, the container surfaces are selected as the interface surfaces, resulting in spherical, cylindrical, rectangular, conic, and disc surface sources. 2.2 Doing the Monte Carlo simulation 2 The complex sources and geometric region is simulated using Monte Carlo method first. The general purpose Monte Carlo code FLUKA [ 17 17 , 18 18 ] is run to get the particle information on the interface. "USERDUMP" is added in the input file to activate the user routine, "MGDRAW", to score the boundary crossing events. The interface is meshed into differential surface elements (Δ A ). The double differential (energy and angle) particle distribution across each surface element is recorded. 2.3 Setting up the interface program 2 The particle information on the interface obtained from the FLUKA simulation should be converted into the surface source of the point kernel calculation. The connection could be established by solving the Boltzmann equation. It could be strictly proved from the Boltzmann equation that the arbitrary source and materials distributed inside the interface could be replaced by the flow rate cross the surface [ 2 2 ]. The response of the detector is the same in both cases. The flow rate passing through the connection surface is equivalent to the intensity of the surface source. <math display="block" id="M1"><mrow><msubsup><mi>j</mi><mtext>n</mtext><mo>+</mo></msubsup><mo stretchy="false">(</mo><msub><mover accent="true"><mi>r</mi><mo>→</mo></mover><mtext>s</mtext></msub><mn>,</mn><mi>E</mi><mn>,</mn><mover accent="true"><mi>Ω</mi><mo>→</mo></mover><mo stretchy="false">)</mo><mo>−</mo><msubsup><mi>j</mi><mtext>n</mtext><mo>−</mo></msubsup><mo stretchy="false">(</mo><msub><mover accent="true"><mi>r</mi><mo>→</mo></mover><mtext>s</mtext></msub><mn>,</mn><mi>E</mi><mn>,</mn><mover accent="true"><mi>Ω</mi><mo>→</mo></mover><mo stretchy="false">)</mo><mo>=</mo><msub><mi>S</mi><mtext>A</mtext></msub><mo stretchy="false">(</mo><msub><mover accent="true"><mi>r</mi><mo>→</mo></mover><mtext>s</mtext></msub><mn>,</mn><mi>E</mi><mn>,</mn><mover accent="true"><mi>Ω</mi><mo>→</mo></mover><mo stretchy="false">)</mo><mn>.</mn></mrow></math> In which S A is the source intensity at the interface, and <math id="M2"><mrow><msubsup><mi>j</mi><mtext>n</mtext><mo>±</mo></msubsup><mo stretchy="false">(</mo><msub><mover accent="true"><mi>r</mi><mo>→</mo></mover><mtext>s</mtext></msub><mn>,</mn><mi>E</mi><mn>,</mn><mover accent="true"><mi>Ω</mi><mo>→</mo></mover><mo stretchy="false">)</mo><mo>=</mo><mover accent="true"><mi>n</mi><mo>^</mo></mover><mo>⋅</mo><mover accent="true"><mi>Ω</mi><mo>→</mo></mover><msup><mi>ϕ</mi><mo>±</mo></msup><mo stretchy="false">(</mo><msub><mover accent="true"><mi>r</mi><mo>→</mo></mover><mtext>s</mtext></msub><mn>,</mn><mi>E</mi><mn>,</mn><mover accent="true"><mi>Ω</mi><mo>→</mo></mover><mo stretchy="false">)</mo></mrow></math> . Where "+" denotes the direction from the Monte Carlo region to the point kernel region, "-" denotes the other way round, <math id="M3"><mover accent="true"><mi>n</mi><mo>^</mo></mover></math> is the outward normal of the surface at any given point, <math id="M4"><mrow><msub><mover accent="true"><mi>r</mi><mo>→</mo></mover><mtext>s</mtext></msub></mrow></math> , E is the particle energy, <math id="M5"><mover accent="true"><mi>Ω</mi><mo>→</mo></mover></math> is the direction vector, and <math id="M6"><mrow><msup><mi>ϕ</mi><mo>±</mo></msup><mo stretchy="false">(</mo><msub><mover accent="true"><mi>r</mi><mo>→</mo></mover><mtext>s</mtext></msub><mn>,</mn><mi>E</mi><mn>,</mn><mover accent="true"><mi>Ω</mi><mo>→</mo></mover><mo stretchy="false">)</mo></mrow></math> denotes the fluence rate. The flow rate from the point kernel region to the Monte Carlo region could be ignored for two reasons: 1) There is no source in the point kernel region. 2) The contributions of the back scattered particles from the reflection surface decreases inversely with the square of the distance. Generally, the interface is far from the reflection surface. Therefore, for the equivalent surface source, <math display="block" id="M7"><mrow><msubsup><mi>j</mi><mtext>n</mtext><mo>−</mo></msubsup><mo stretchy="false">(</mo><msub><mover accent="true"><mi>r</mi><mo>→</mo></mover><mtext>s</mtext></msub><mn>,</mn><mi>E</mi><mn>,</mn><mover accent="true"><mi>Ω</mi><mo>→</mo></mover><mo stretchy="false">)</mo><mo>=</mo><mn>0,</mn><msubsup><mi>j</mi><mtext>n</mtext><mo>+</mo></msubsup><mo stretchy="false">(</mo><msub><mover accent="true"><mi>r</mi><mo>→</mo></mover><mtext>s</mtext></msub><mn>,</mn><mi>E</mi><mn>,</mn><mover accent="true"><mi>Ω</mi><mo>→</mo></mover><mo stretchy="false">)</mo><mo>=</mo><msub><mi>S</mi><mtext>A</mtext></msub><mo stretchy="false">(</mo><msub><mover accent="true"><mi>r</mi><mo>→</mo></mover><mtext>s</mtext></msub><mn>,</mn><mi>E</mi><mn>,</mn><mover accent="true"><mi>Ω</mi><mo>→</mo></mover><mo stretchy="false">)</mo><mn>.</mn></mrow></math> The angular distribution of the surface source is usually not isotropic, but we can use a more generalized form to express it as a power distribution of cos θ [ 2 2 ], <math display="block" id="M8"><mrow><mi>f</mi><mo stretchy="false">(</mo><mi>θ</mi><mo stretchy="false">)</mo><mo>=</mo><mfrac><mn>1</mn><mrow><mn>2</mn><mi>π</mi></mrow></mfrac><msubsup><mi>j</mi><mtext>n</mtext><mo>+</mo></msubsup><mo stretchy="false">(</mo><msub><mover accent="true"><mi>r</mi><mo>→</mo></mover><mtext>s</mtext></msub><mn>,</mn><mi>E</mi><mn>,</mn><mover accent="true"><mi>Ω</mi><mo>→</mo></mover><mo stretchy="false">)</mo><mo stretchy="false">(</mo><mi>m</mi><mo>+</mo><mn>1</mn><mo stretchy="false">)</mo><msup><mstyle mathsize="140%" displaystyle="true"><mrow><mi>cos</mi></mrow></mstyle><mi>m</mi></msup><mi>θ</mi><mn>,</mn></mrow></math> where θ is the angle between the emission and the outward normal of the surface element, and m is a variable parameter, which can be found using the least-squares fitting. The integration of the function f ( θ ) over solid angle is equal to <math id="M9"><mrow><msubsup><mi>j</mi><mtext>n</mtext><mo>+</mo></msubsup><mo stretchy="false">(</mo><msub><mover accent="true"><mi>r</mi><mo>→</mo></mover><mtext>s</mtext></msub><mn>,</mn><mi>E</mi><mn>,</mn><mover accent="true"><mi>Ω</mi><mo>→</mo></mover><mo stretchy="false">)</mo></mrow></math> , <math display="block" id="M10"><mrow><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mi>Ω</mi></msub></mrow></mstyle><mi>f</mi><mo stretchy="false">(</mo><mi>θ</mi><mo stretchy="false">)</mo><mtext>d</mtext><mi>Ω</mi><mo>=</mo><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mn>0</mn><mrow><mn>2</mn><mi>π</mi></mrow></msubsup></mrow></mstyle><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mn>0</mn><mrow><mi>π</mi><mo>/</mo><mn>2</mn></mrow></msubsup></mrow></mstyle><mi>f</mi><mo stretchy="false">(</mo><mi>θ</mi><mo stretchy="false">)</mo><mi>sin</mi><mi>θ</mi><mtext>d</mtext><mi>θ</mi><mtext>d</mtext><mi>ϕ</mi><mo>=</mo><msubsup><mi>j</mi><mtext>n</mtext><mo>+</mo></msubsup><mo stretchy="false">(</mo><msub><mover accent="true"><mi>r</mi><mo>→</mo></mover><mtext>s</mtext></msub><mn>,</mn><mi>E</mi><mn>,</mn><mover accent="true"><mi>Ω</mi><mo>→</mo></mover><mo stretchy="false">)</mo><mn>,</mn></mrow></math> where the distribution is supposed to be symmetric with the azimuthal angle, ϕ . For an arbitrarily source, the angular distribution of the emissions would not be expressed as a continuous function of cos θ . The discrete method is used instead in this situation. 2.4 Doing the point kernel approach 2 After getting the angular distribution on the interface, the point kernel approach is done. The steps to calculate the detector response are summarized as follows: • Use the same mesh of the interface as Sec. 2.2, • Calculate the detector response to the surface elements, • Add the contribution of all the surface elements. The most important part of doing the point kernel approach is to calculate the detector response to the surface elements. With a plane surface source as an example, Fig. 1 Fig. 1 shows the geometry for the calculation. For an arbitrary surface element, dS , at point A , the detector is located at point P . The projection of P on the plane is P’. The angle between <math id="M11"><mrow><mover accent="true"><mi>A</mi><mo>→</mo></mover><mi>P</mi></mrow></math> and the outward normal of the surface element around point A , is characterized by θ . The azimuthal angle of detector P is expressed as ϕ . θ and ϕ could be calculated using <math display="block" id="M12"><mrow><mi>cos</mi><mi>θ</mi><mo>=</mo><mfrac><mrow><mover accent="true"><mi>A</mi><mo>→</mo></mover><mi>P</mi><mo>⋅</mo><mover accent="true"><mi>z</mi><mo>^</mo></mover></mrow><mrow><mo>|</mo><mover accent="true"><mi>A</mi><mo>→</mo></mover><mi>P</mi><mo>|</mo></mrow></mfrac><mn>,</mn><mi>cos</mi><mi>ϕ</mi><mo>=</mo><mfrac><mrow><mover accent="true"><mi>A</mi><mo>→</mo></mover><msup><mi>P</mi><mo>′</mo></msup><mo>⋅</mo><mover accent="true"><mi>x</mi><mo>^</mo></mover></mrow><mrow><mo>|</mo><mover accent="true"><mi>A</mi><mo>→</mo></mover><msup><mi>P</mi><mo>′</mo></msup><mo>|</mo></mrow></mfrac><mn>.</mn></mrow></math> For the case when the emitting angle is symmetric with ϕ , it could be expressed as a cosine function of θ . The detector response is expresses as: <math display="block" id="M13"><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mi>d</mi><mi>R</mi><mo>=</mo></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mfrac><mn>1</mn><mrow><mn>2</mn><mi>π</mi></mrow></mfrac><msup><mi>r</mi><mrow><mo>−</mo><mn>2</mn></mrow></msup><mi>B</mi><mo stretchy="false">(</mo><mi>E</mi><mn>,</mn><mi>d</mi><mo stretchy="false">)</mo><msup><mi>e</mi><mrow><mo>−</mo><mi>μ</mi><mi>d</mi><mo stretchy="false">(</mo><mi>θ</mi><mo stretchy="false">)</mo></mrow></msup><mi>ℜ</mi><mo stretchy="false">(</mo><mi>E</mi><mo stretchy="false">)</mo><msubsup><mi>j</mi><mi>n</mi><mo>+</mo></msubsup><mo stretchy="false">(</mo><mi>m</mi><mo>+</mo><mn>1</mn><mo stretchy="false">)</mo><msup><mstyle mathsize="140%" displaystyle="true"><mrow><mi>cos</mi></mrow></mstyle><mi>m</mi></msup><mi>θ</mi><mi>d</mi><mi>E</mi><mi>d</mi><mi>S</mi><mn>.</mn></mrow></mtd></mtr></mtable></mrow></math> where <math id="M14"><mrow><mi>r</mi><mo>=</mo><mo>|</mo><mover accent="true"><mi>A</mi><mo>→</mo></mover><mi>P</mi><mo>|</mo></mrow></math> , B ( E , d ) is the gamma-ray buildup factor through the shield, d ( θ ) is the shield thickness along the θ direction, μ is the photon attenuation coefficient, and R ( E ) is the detector response function. The widely used ANSI/ANS-6.4.3-1991 buildup factors and attenuation coefficients are adopted [ 19 19 ] in the paper. Since the buildup factors are provided for point isotropic sources in a homogeneous infinite medium, the error increases for practical situations, e.g. a shield of finite extent and slant-penetration problems. Situation-specific buildup factors or situation-specific corrections should be used for the special cases. For example, different series of a buildup factor could be used for a parallel beam obliquely incident on a slab [ 20 20 ]. For an arbitrary emitting angle, which is depending not only on θ but also on ϕ , the discrete method is used. The detector response is: <math display="block" id="M15"><mrow><mi>d</mi><mi>R</mi><mo>=</mo><msup><mi>r</mi><mrow><mo>−</mo><mn>2</mn></mrow></msup><msub><mi>H</mi><mrow><msub><mi>θ</mi><mi>i</mi></msub><mn>,</mn><msub><mi>ϕ</mi><mi>j</mi></msub></mrow></msub><mi>B</mi><mo stretchy="false">(</mo><mi>E</mi><mn>,</mn><mi>d</mi><mo stretchy="false">)</mo><msup><mi>e</mi><mrow><mo>−</mo><mi>μ</mi><mi>d</mi><mo stretchy="false">(</mo><msub><mi>θ</mi><mi>i</mi></msub><mn>,</mn><msub><mi>ϕ</mi><mi>j</mi></msub><mo stretchy="false">)</mo></mrow></msup><mi>d</mi><mi>E</mi><mi>d</mi><mi>S</mi><mn>.</mn></mrow></math> where H θ i , ϕ j is the dose equivalent at one unit distance from the source with no shield along ( θ i , ϕ j ) direction. Integrating the differential detector response over energy and the whole interface, the total detector response is obtained. <math display="block" id="M16"><mrow><msub><mi>R</mi><mrow><mtext>total</mtext></mrow></msub><mo>=</mo><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mi>S</mi></msub></mrow></mstyle><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mrow><mi>E</mi><mo>=</mo><mn>0</mn></mrow><mrow><msub><mi>E</mi><mrow><mtext>max</mtext></mrow></msub></mrow></msubsup></mrow></mstyle><mtext>d</mtext><mi>R</mi><mn>.</mn></mrow></math> 2.5 Considering the wall/floor reflection 2 The reflection from the floor and the walls are considered by introducing a reflection factor, which can be obtained from the albedo concept. The seven-parameter approximation of the dose albedo given by Chilton and Huddleston [ 2 2 ] is expressed as <math display="block" id="M17"><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><msub><mi>α</mi><mtext>D</mtext></msub><mo stretchy="false">(</mo><msub><mi>E</mi><mn>0</mn></msub><mn>,</mn><msub><mi>θ</mi><mn>0</mn></msub><mo>;</mo><mi>θ</mi><mn>,</mn><mi>ϕ</mi><mo stretchy="false">)</mo></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mo>=</mo><mi>F</mi><mo stretchy="false">(</mo><msub><mi>E</mi><mn>0</mn></msub><mn>,</mn><msub><mi>θ</mi><mn>0</mn></msub><mo>;</mo><mi>θ</mi><mn>,</mn><mi>ϕ</mi><mo stretchy="false">)</mo><mfrac><mrow><mi>C</mi><mo stretchy="false">(</mo><msub><mi>E</mi><mn>0</mn></msub><mo stretchy="false">)</mo><mi>K</mi><mo stretchy="false">(</mo><msub><mi>E</mi><mn>0</mn></msub><mn>,</mn><msub><mi>θ</mi><mi>s</mi></msub><msup><mrow><mo stretchy="false">)</mo><mn>10</mn></mrow><mrow><mn>26</mn></mrow></msup><mo>+</mo><msup><mi>C</mi><mo>'</mo></msup><mo stretchy="false">(</mo><msub><mi>E</mi><mn>0</mn></msub><mo stretchy="false">)</mo></mrow><mrow><mn>1</mn><mo>+</mo><mo stretchy="false">(</mo><mi>cos</mi><msub><mi>θ</mi><mn>0</mn></msub><mo>/</mo><mi>cos</mi><mi>θ</mi><mo stretchy="false">)</mo><mo stretchy="false">(</mo><mn>1</mn><mo>+</mo><mn>2</mn><msub><mi>E</mi><mn>0</mn></msub><mtext>vers</mtext><msub><mi>θ</mi><mi>s</mi></msub><msup><mo stretchy="false">)</mo><mrow><mn>1</mn><mo>/</mo><mn>2</mn></mrow></msup></mrow></mfrac><mn>,</mn></mrow></mtd></mtr></mtable></mrow></math> in which the factor F is given by <math display="block" id="M18"><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mi>F</mi><mo stretchy="false">(</mo><msub><mi>E</mi><mn>0</mn></msub><mn>,</mn><msub><mi>θ</mi><mn>0</mn></msub><mo>;</mo><mi>θ</mi><mn>,</mn><mi>ϕ</mi><mo stretchy="false">)</mo><mo>=</mo><msub><mi>A</mi><mn>1</mn></msub><mo stretchy="false">(</mo><msub><mi>E</mi><mn>0</mn></msub><mo stretchy="false">)</mo><mo>+</mo><msub><mi>A</mi><mn>2</mn></msub><mo stretchy="false">(</mo><msub><mi>E</mi><mn>0</mn></msub><mo stretchy="false">)</mo><msup><mrow><mtext>vers</mtext></mrow><mn>2</mn></msup><msub><mi>θ</mi><mn>0</mn></msub><mo>+</mo><msub><mi>A</mi><mn>3</mn></msub><mo stretchy="false">(</mo><msub><mi>E</mi><mn>0</mn></msub><mo stretchy="false">)</mo><msup><mrow><mtext>vers</mtext></mrow><mn>2</mn></msup><mi>θ</mi></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mo>+</mo><msub><mi>A</mi><mn>4</mn></msub><mo stretchy="false">(</mo><msub><mi>E</mi><mn>0</mn></msub><mo stretchy="false">)</mo><msup><mrow><mtext>vers</mtext></mrow><mn>2</mn></msup><msub><mi>θ</mi><mn>0</mn></msub><msup><mrow><mtext>vers</mtext></mrow><mn>2</mn></msup><mi>θ</mi><mo>+</mo><msub><mi>A</mi><mn>5</mn></msub><mo stretchy="false">(</mo><msub><mi>E</mi><mn>0</mn></msub><mo stretchy="false">)</mo><mtext>vers</mtext><msub><mi>θ</mi><mn>0</mn></msub><mtext>vers</mtext><mi>θ</mi><mtext>vers</mtext><mi>ϕ</mi><mn>,</mn></mrow></mtd></mtr></mtable></mrow></math> where E 0 is the energy of the incident particle, θ 0 and θ are the incident and emergent polar angle, respectively, with respect to the outward normal of the reflecting surface, vers θ =1-cos θ , ϕ is the emergent azimuthal angle, C , C’, A 1 , A 2 , A 3 , A 4 , and A 5 are fitted parameters, which depend not only on energy, but also on the composition of the reflecting medium [ 21 21 ], and K ( E 0 , θ s ) is the Klein-Nishina energy scattering cross section for scattering angle θ s in the units of square centimeters. The response of the detector, which is located at the direction ( θ , ϕ ) with respect to the reflecting surface dS’ , is <math display="block" id="M19"><mrow><mi>d</mi><msub><mi>R</mi><mi>r</mi></msub><mo>=</mo><msub><mi>R</mi><mn>0</mn></msub><mfrac><mrow><msub><mi>α</mi><mtext>D</mtext></msub><mo stretchy="false">(</mo><msub><mi>E</mi><mn>0</mn></msub><mn>,</mn><msub><mi>θ</mi><mn>0</mn></msub><mo>;</mo><mi>θ</mi><mn>,</mn><mi>ϕ</mi><mo stretchy="false">)</mo><mi>cos</mi><msub><mi>θ</mi><mn>0</mn></msub><mi>d</mi><msup><mi>S</mi><mo>′</mo></msup></mrow><mrow><msubsup><mi>r</mi><mn>2</mn><mn>2</mn></msubsup></mrow></mfrac><mn>,</mn></mrow></math> where R 0 represents the response at the location of dS’ , due to incident particles, and r 2 is the distance from the detector to dS’ . Integrating over the area of the reflecting surface, the total reflected dose is obtained. It should be noticed that as the differential surface dS’ changes, the variables θ 0 , θ , ϕ , and r 2 change as well. 2.6 Getting the 3D radiation field 2 Figure 2 Figure 2 shows the flow chart of the coupled model. The interface program and the point kernel approach module are developed using the Fortran language. Combing the Monte Carlo and point kernel results, the 3D radiation field is obtained. 3 Verification of the coupled program 1 Two models from simplified typical waste storage facilities are used to verify the coupled program. One is used to test the case when the emitting angle from the equivalent surface source is a power distribution of cos θ , and the other one is designed for an arbitrary emitting angle, which is depending not only on θ but also on ϕ . 3.1 Case A: the emitting angle is a power distribution of cos θ 2 3.1.1 The geometry and source 3 As depicted in Fig. 3(a) Fig. 3(a) , there are two rooms in the facility. Two identical vessels containing 1.5 g/cm 3 silt with radioactive nuclide Cs-137 are placed in the room on the right side. The walls of the vessels are composed by 10 mm steel. A concrete wall of 20 cm at x =-673●thin;cm separates the vessels from the room on the left side. The C-C’ and D-D’ cross sections of the vessels, which are identical, are shown by Fig. 3(b) Fig. 3(b) . The major part of the vessel is a cylinder that is 3.1 meters high with an inverted frustum of cone of 1.0 meter on the bottom. The source is exponentially distributed throughout the z direction and inclined to the negative direction along the x -axis. 3.1.2 The coupled computation 3 A full Monte Carlo simulation is conducted first using FLUKA to do the benchmark. The "BIASING" card is used to compensate the photon attenuation caused by the concrete wall and by the geometric distance from the source. The surfaces of the vessels are selected as the interface for the coupled computation, resulting in cylindrical, conic, and disc surface sources. The cylindrical and conic surface sources are meshed in cylindrical coordinates, while the disc surface source is meshed in Cartesian coordinates. The mesh size for each element is selected to make sure the deviation of the final results is less than 2% when increasing the mesh number. The entire energy range is subdivided into 15 groups and the angular variable is discretized into 15 parts. To obtain the equivalent surface source, the current <math id="M20"><mrow><msubsup><mi>j</mi><mtext>n</mtext><mo>+</mo></msubsup></mrow></math> on each surface segment is given by the FLUKA simulation. Figure 4 Figure 4 shows the double differential spectra of emissions with θ in different energy intervals for a selected frustum surface element. Three energy intervals are taken for drawing comparison, where E 1 ∈[0.596,0.662] MeV, E 2 ∈[0.529,0.596] MeV, E 3 ∈[0.463,0.529] MeV. The distribution of the emissions is expressed as the power of cos θ , <math id="M21"><mrow><msubsup><mi>j</mi><mtext>n</mtext><mo>+</mo></msubsup><mo stretchy="false">(</mo><mi>m</mi><mo>+</mo><mn>1</mn><mo stretchy="false">)</mo><msup><mstyle mathsize="140%" displaystyle="true"><mrow><mi>cos</mi></mrow></mstyle><mi>m</mi></msup><mi>θ</mi><mo>/</mo><mn>2</mn><mi>π</mi></mrow></math> . The parameter, m , is equal to 2, which is obtained using the least-squares fitting. It is shown by the figure that the data fits well with the formula. The point kernel approach is done following the steps of Sec. 2.4, using the current <math id="M22"><mrow><msubsup><mi>j</mi><mtext>n</mtext><mo>+</mo></msubsup></mrow></math> given by the FLUKA simulation. The reflections caused by the wall near the vessels and the floor below them are considered. The incident angle, θ 0 , is selected to be zero due to the approximate normal incidence. 3.1.3 The verification of the coupled program 3 To verify the new program, the results of coupled computation are compared with the full FLUKA simulation. The ambient dose, equivalent as a function of x is shown in Fig. 5 Fig. 5 for the line of z =250●thin;cm, y =0●thin;cm and z =250●thin;cm, y =235●thin;cm, where y =0●thin;cm is the symmetrical plane of geometry, and y =235●thin;cm is the plane passing through the central axis of a vessel as depicted in Fig. 3 Fig. 3 . The solid black line shows the large divergence of the point kernel approach at the position which has the label of "Container Surface". It is an intrinsic drawback of the point kernel method. When the detector is very "close" to the source, the term 1/ r 2 in Eq. 6 and Eq. 7 approaches infinity. In such cases, the results oscillate substantially when the mesh of the interface is changed slightly. The radiation field in such region is selected from the FLUKA simulation in the coupled method. The dose equivalent decreases fast when passing through the concrete wall located from x =-663●thin;cm to x =-683●thin;cm (called "inwall" for short). The selection of the buildup factor is always an artistic problem in the practical design, especially for odd shaped or large objects. The buildup factor for an isotropic point source in an infinite homogeneous medium is usually chosen in the shielding calculations, since it is likely to be a conservative estimation. However, as shown in Fig. 4 Fig. 4 , since most of the particles move forward along the outward normal of the surface, the beam approximately parallel illuminates the "inwall" in our case. Therefore, the buildup factor for the parallel beam fits better with the FLUKA simulation. 3.2 Case B: the emitting angle is a discrete distribution 2 3.2.1 The geometry and source 3 As depicted in Fig. 6 Fig. 6 , three identical vessels the same as the one in Sec. 3.1 are placed on the right side of the room. There are two 20 cm concrete walls named by "wall1" and "wall2". The dotted line box enclosing the vessels is selected as the interface. 3.2.2 The coupled computation 3 The steps to do the coupled computation are similar to Sec. 3.1. The biggest difference is: since the emitting angle from the interface is depending not only on θ but also on ϕ , the discrete method is used. For a given detector, the response from each side of the interface is added. The buildup factor for the parallel beam is chosen in the calculation. Figure 7 Figure 7 shows the ambient dose equivalent as a function of x , with z =250●thin;cm. Two representative lines, y =-234●thin;cm and y =403●thin;cm, are selected to compare with the FLUKA simulation. The two results match very well on a large scale for the ambient dose equivalent, which drops about five orders of magnitude. We also see large fluctuation for the FLUKA results after passing "wall2", although 1.2× 10 10 particles have been simulated. 3.3 The advantage of the coupled scheme 2 The good agreement of the results in Sec. 3.1 and Sec. 3.2 indicate the reliability of the coupled Monte Carlo-Point Kernel scheme. Comparing with the Monte Carlo method, the coupled scheme can speed up the calculation under the premise of ensuring the accuracy. For case A, using the desktop computer with Intel(R) Core(TM) i3-3240 CPU @3.4GHz, 4GB memory, for a mesh of 70× 50× 20, it takes eight hours for the full FLUKA simulation to get the relative uncertainties of 5% in the right corner and 15% in the left corner of the left room. The CPU time to get the equivalent surface source in the coupled technique is twenty minutes within the relative uncertainty of 1%, and the time for the point kernel calculation is thirty-five minutes using the same mesh size. For case B, using the same computer and the same mesh, thirty hours were spent for the full FLUKA simulation due to the two thick shields, although the "BIASING" card is used to reduce the variance. The relative uncertainty is 12% in the right corner and 24% in the left corner of the most left room for the FLUKA simulation. The coupled technique takes seventy-five minutes to get the equivalent surface source with the relative uncertainty of 5% and two hours and fourteen minutes to do the point kernel calculation. The results show that the coupled technique saves approximately one order of time compared to the full Monte Carlo simulation. 4 Conclusion 1 The coupled Monte Carlo-Point Kernel scheme, which could meet the requirements of the fast radiation field calculation for large-scale nuclear installations, has been presented and validated. The coupled program system takes the advantages of both the Monte Carlo simulation and the point kernel method. The combination of these features makes it possible to deal with the complex geometry and source using the Monte Carlo model and do fast dose assessments by the point kernel approach in the large-scale bulk shielding region. Two models from simplified typical waste storage facilities are selected to do the code verification. Although the test models are rather simple, they have the following features: • The source is non-standard; • The shield is bulk; • The building house has large dimension; • A fast 3D dose field map is needed. As the large-scale complex nuclear installations share the same features as pointed above, it can be concluded that the simple models test and verify the ability of the coupled scheme to do such simulations. The coupled scheme will be extended to be capable of dealing with more complex cases in the follow up study. Using the actual scene of the large-scale nuclear facilities, the source and the geometry will inevitably be more complex. Several examples covering a wide range of expected configurations would be examined to determine the overall validity of the method for use during the decommissioning process. References J. B. Thevenon , NARVEOS: a New Tool Supporting ALARA Studies . IAEA Training Course on Decommissioning Dose Assessment & Dose Optimization, MARCOULE CEA Center , ( 2011 ) NARVEOS: a New Tool Supporting ALARA Studies B. Chilton , J. K. Shultis , R. E. Faw , Principles of Radiation Shielding . Englewood Cliffs, New Jersey 07632: Prentice-Hall, Inc. , 161-169, 266-284, 310-365 ( 1984 ) Principles of Radiation Shielding NATIONAL COUNCIL ON RADIATION PROTECTION AND MEASUREMENTS , Radiation Protection for Particle Accelerator Facilities. NCRP Report No. 144 , ( 2003 ) Radiation Protection for Particle Accelerator Facilities. NCRP Report No. 144 J. C. Wagner , D. E. Peplow , S. W. Mosher , FW-CADIS Method For Global and Regional Variance Reduction Of Monte Carlo Radiation Transport Calculations . Nucl. Sci. Eng . 176 , 37 - 57 ( 2014 ). doi: 10.13182/NSE12-33 http://doi.org/10.13182/NSE12-33 FW-CADIS Method For Global and Regional Variance Reduction Of Monte Carlo Radiation Transport Calculations J. C. Wagner , A. Haghighat , Automated Variance Reduction of Monte Carlo Shielding Calculations Using the Discrete Ordinates Adjoint Function . Nucl. Sci. Eng . 128 : 186 - 208 ( 1998 ). doi: 10.13182/NSE98-2 http://doi.org/10.13182/NSE98-2 Automated Variance Reduction of Monte Carlo Shielding Calculations Using the Discrete Ordinates Adjoint Function F. Vermeersch , Dose Assessment and Dose Optimization in Decommissioning Using the VISIPLAN 3D ALARA Planning Tool . Radiation Protection and Decommissioning Conference ABR/BVS , Brussels , ( 2003 ) Dose Assessment and Dose Optimization in Decommissioning Using the VISIPLAN 3D ALARA Planning Tool I. Szőke , M. N. Louka , T. R. Bryntesen , et al ., Real-time 3D radiation risk assessment supporting simulation of work in nuclear environments . J. Radiol. Prot . 34 ( 2 ), 389 - 416 ( 2014 ). doi: 10.1088/0952-4746/34/2/389 http://doi.org/10.1088/0952-4746/34/2/389 Real-time 3D radiation risk assessment supporting simulation of work in nuclear environments Electric Power Research Institute (EPRI), Work Planning and 3D Visualization Tool . www.epri.com www.epri.com Electric Power Research Institute (EPRI), Work Planning and 3D Visualization Tool P. Tran , S. Bickford , Collaborative Project: 3D EDE Radiation Exposure Planning Tool . Plant Productivity Workshop , ( 2010 ) Collaborative Project: 3D EDE Radiation Exposure Planning Tool Siemens Industry Software , ALARA Planning: Siemens PLM Software and Microsoft Engineer a Partnership to Ensure Employee Health and Safety . www.siemens.com/plm/alara www.siemens.com/plm/alara ALARA Planning: Siemens PLM Software and Microsoft Engineer a Partnership to Ensure Employee Health and Safety F. X. Gallmeier , R. E. Pevey , Creation of a Set of Interface Utilities to Allow Coupled Monte Carlo/Discrete Ordinate Shielding Analysis . Third International Topical Meeting on Nuclear Applications of Accelerator Technology , USA : American Nuclear Society , ( 1999 ) Creation of a Set of Interface Utilities to Allow Coupled Monte Carlo/Discrete Ordinate Shielding Analysis J. O. Johnson , R. T. Santoro , R. A. Lillie , et al ., The SNS Target Station Preliminary Title I: Shielding Analyses . J. Nucl. Sci. Technol . 37 ( sup1 ), 35 - 39 ( 2000 ). doi: 10.1080/00223131.2000.10874842 http://doi.org/10.1080/00223131.2000.10874842 The SNS Target Station Preliminary Title I: Shielding Analyses Y. Chen , Coupled Monte Carlo Discrete Ordinates Computational Scheme for Three-Dimensional Shielding Calculations of Large and Complex Nuclear Facilities. Ph.D Thesis , Forschungszentrum Karlsruhe GmbH , karlsruhe, Germany , ( 2005 ) Coupled Monte Carlo Discrete Ordinates Computational Scheme for Three-Dimensional Shielding Calculations of Large and Complex Nuclear Facilities. Ph.D Thesis J. R. Han , Q. F. Liu , H. Y. Chen , et al ., Development of three dimensional discrete ordinates-Monte Carlo coupled system . Nucl. Sci. Tech . 25 ( S1 ), S010601 ( 2014 ) doi: 10.13538/j.1001-8042/nst.25.S010601 http://doi.org/10.13538/j.1001-8042/nst.25.S010601 Development of three dimensional discrete ordinates-Monte Carlo coupled system J. K. Shultis , R. E. Faw , Radiation Shielding Technology . Health Phys . 88 ( 6 ), 587 - 612 ( 2005 ). doi: 10.1097/01.HP.0000148615.73825.b1 http://doi.org/10.1097/01.HP.0000148615.73825.b1 Radiation Shielding Technology Y. J. Bi , H. W. Yang , J. S. Guo , et al ., Study on Fast Radiation Transport Calculation Methodology in Decommissioning of Research Reactor . Atomic Energy Science and Technology , 49 ( sup1 ): 470 - 474 ( 2015 ). (in Chinese) Study on Fast Radiation Transport Calculation Methodology in Decommissioning of Research Reactor G. Battistoni , S. Muraro , P.R. Sala , et al ., The FLUKA code: Description and benchmarking, Proceedings of the Hadronic Shower Simulation Workshop. Fermilab, (2006) . M. Albrow , R. Raja eds., AIP Conference Proceeding 896 , 31 - 49 ( 2007 ) The FLUKA code: Description and benchmarking, Proceedings of the Hadronic Shower Simulation Workshop. Fermilab, (2006) A. Fasso , A. Ferrari , J. Ranft , et al ., FLUKA: a Multi-Particle Transport Code. CERN-2005-10, INFN/TC_05/11, SLAC-R-773 , ( 2005 ) FLUKA: a Multi-Particle Transport Code. CERN-2005-10, INFN/TC_05/11, SLAC-R-773 ANSI/ANS-6.4.3-1991, Gamma-Ray Attenuation Coefficients and Buildup Factors for Engineering Materials . La Grange Park : American Nuclear Society , ( 1991 ) E. M. Fournie , A. B. Chilton , Gamma-ray Buildup Factors for Concrete Slab Shields under Slant Incidence Conditions . Nucl. Sci. Eng . 76 : 66 - 69 ( 1980 ). doi: 10.13182/NSE80-2 http://doi.org/10.13182/NSE80-2 Gamma-ray Buildup Factors for Concrete Slab Shields under Slant Incidence Conditions R. C. Brockhoff , J. K. Shultis , Data for the Calculation of Albedos from Concrete, Iron, Lead and Water for Photons and Neutrons . Kansas State University, Engineering Experiment Station , Report 301 , ( 2005 ) Data for the Calculation of Albedos from Concrete, Iron, Lead and Water for Photons and Neutrons 10.1007/s41365-017-0219-6 This work was Supported by the National Natural Science Foundation of China (No.11405277) 2016-11-04 2016-12-16 2016-12-25 2017-03-29 2017-05-01 2021-04-11 © China Science Publishing & Media Ltd. (Science Press), Shanghai Institute of Applied Physics, the Chinese Academy of Sciences, ChineseNuclear Society and Springer Nature Singapore Pte Ltd. 2017 This is a post-peer-review, pre-copyedit version of an article published in Nuclear Science and Techniques. The final authenticated version is available online at: http://dx.doi.org/10.1007/s41365-017-0219-6 http://dx.doi.org/10.1007/s41365-017-0219-6 . 2017 Nuclear Science and Techniques 05 28