Development of three dimensional discrete ordinates-Monte Carlo coupled system

Development of three dimensional discrete ordinates-Monte Carlo coupled system 1 corresp first-author Corresponding author, hanjingru@163.com Corresponding author, hanjingru@163.com hanjingru@163.com 1 1 1 Nuclear and Radiation Safety Center of Ministry of Environmental Protection, Beijing100082, China Nuclear and Radiation Safety Center of Ministry of Environmental Protection, Beijing100082, China The radiation shielding calculations of large nuclear facilities are complicated due to their bulk shields and complex geometries. A program system based on three-dimensional discrete ordinates ( S N )-Monte Carlo (MC) coupled method has been developed to solve this kind of shielding problems. S N method is used to treat the deep penetration problem in the bulk shield, and MC simulation is performed with complex geometry. The coupled method is implemented in the interface program to calculate the probability of particle variables from the S N angular flux distribution. A special source routine is written and linked to MC code to processed these probability arrays and sample the particle variables to generate source for the use of MC code. The H. B. Robinson Unit 2 pressure vessel benchmark was used to validate the S N -MC coupled program system. The specific activities for the six dosimeters in the surveillance capsule were calculated. Satisfactory agreements were obtained from the comparisons of S N -MC results with those obtained from measurement and other computer codes. The result demonstrates that the S N -MC coupling scheme with the program system is suitable to treat three dimensional shielding problems with satisfactory accuracy. Discrete ordinates Monte Carlo Coupled Benchmark validation Shielding I. INTRODUCTION 1 Radiation shielding calculation of large scale nuclear facilities is complicated due to the bulk shield and complex geometry. Discrete ordinates ( S N ) and Monte Carlo (MC) methods are commonly used in the design of radiation shielding. MC method is able to treat the complex geometry, but time-consuming in dealing with the deep penetration problems. S N method has great computational efficiency, but it is quite costly in computer memory. Therefore, for the shielding analysis of a large and complex system, the coupled method which combines the advantages of both the MC and S N approaches is very suitable. The MC method is used to simulate the particle transport in the complex geometry, and the S N method is used to treat the deep penetration problems in the bulk shield. The existing S N -MC coupled codes include DOT-DOMINO-MORSE [ 1 1 ] and DORT-PROBGEN-MCNP [ 2 2 ]. However, these codes cannot support full three-dimensional geometry. Masahiko Kurosawa developed TORT-MCNP [ 3 3 ] for the calculation of neutron flux around a core of BWR, but it only for Cartesian geometry. In this work, a program system has been devised to make full three-dimensional S N -MC coupling calculations in both Cartesian and Cylindrical geometry. The validation calculation of the H. B. Robinson Unit 2 (HBR-2) pressure vessel benchmark problem [ 4 4 ] is presented and discussed. II. S N -MC COUPLED METHOD 1 To perform the S N -MC calculation, the model is divided into two parts. The one part of the model is the bulk shield region with simple geometry, and the other part is the complex geometry. The discrete ordinates method has advantage in treating the deep penetration problem with simple geometry. The MC method is able to treat the complex geometry, but time-consuming in dealing with the deep penetration problem. In order to eliminate the weakness of the single method, an interface program is developed to transform the angular fluxes from S N solution into normalized probability distributions. The normalized probability distributions calculated by interface program are as follows: <math display="block" id="M1"><mrow><mi>p</mi><mo stretchy="false">(</mo><mi>i</mi><mo stretchy="false">)</mo><mo>=</mo><mfrac><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>g</mi><mo>−</mo><mn>1</mn></mrow><mi>G</mi></munderover></mstyle><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>m</mi><mo>−</mo><mn>1</mn></mrow><mi>M</mi></munderover></mstyle><msub><mi>w</mi><mi>m</mi></msub><mi>ψ</mi><mo stretchy="false">(</mo><msub><mi>r</mi><mi>i</mi></msub><mn>,</mn><msub><mi>E</mi><mi>g</mi></msub><mn>,</mn><msub><mi>Ω</mi><mi>m</mi></msub><mo stretchy="false">)</mo><mo>|</mo><msub><mi>λ</mi><mi>m</mi></msub><mo>|</mo><msub><mi>A</mi><mi>i</mi></msub></mrow><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow><mi>I</mi></munderover></mstyle><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>g</mi><mo>−</mo><mn>1</mn></mrow><mi>G</mi></munderover></mstyle><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>m</mi><mo>−</mo><mn>1</mn></mrow><mi>M</mi></munderover></mstyle><msub><mi>w</mi><mi>m</mi></msub><mi>ψ</mi><mo stretchy="false">(</mo><msub><mi>r</mi><mi>i</mi></msub><mn>,</mn><msub><mi>E</mi><mi>g</mi></msub><mn>,</mn><msub><mi>Ω</mi><mi>m</mi></msub><mo stretchy="false">)</mo><mo>|</mo><msub><mi>λ</mi><mi>m</mi></msub><mo>|</mo><msub><mi>A</mi><mi>i</mi></msub></mrow></mfrac><mn>,</mn></mrow></math> <math display="block" id="M2"><mrow><mi>p</mi><mo stretchy="false">(</mo><mi>g</mi><mo>/</mo><mi>i</mi><mo stretchy="false">)</mo><mo>=</mo><mfrac><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>m</mi><mo>−</mo><mn>1</mn></mrow><mi>M</mi></munderover></mstyle><msub><mi>w</mi><mi>m</mi></msub><mi>ψ</mi><mo stretchy="false">(</mo><msub><mi>r</mi><mi>i</mi></msub><mn>,</mn><msub><mi>E</mi><mi>g</mi></msub><mn>,</mn><msub><mi>Ω</mi><mi>m</mi></msub><mo stretchy="false">)</mo><mo>|</mo><msub><mi>λ</mi><mi>m</mi></msub><mo>|</mo><msub><mi>A</mi><mi>i</mi></msub></mrow><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>g</mi><mo>−</mo><mn>1</mn></mrow><mi>G</mi></munderover></mstyle><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>m</mi><mo>−</mo><mn>1</mn></mrow><mi>M</mi></munderover></mstyle><msub><mi>w</mi><mi>m</mi></msub><mi>ψ</mi><mo stretchy="false">(</mo><msub><mi>r</mi><mi>i</mi></msub><mn>,</mn><msub><mi>E</mi><mi>g</mi></msub><mn>,</mn><msub><mi>Ω</mi><mi>m</mi></msub><mo stretchy="false">)</mo><mo>|</mo><msub><mi>λ</mi><mi>m</mi></msub><mo>|</mo><msub><mi>A</mi><mi>i</mi></msub></mrow></mfrac><mn>,</mn></mrow></math> <math display="block" id="M3"><mrow><mi>p</mi><mo stretchy="false">(</mo><mi>m</mi><mo>/</mo><mi>g</mi><mo>/</mo><mi>i</mi><mo stretchy="false">)</mo><mo>=</mo><mfrac><mrow><msub><mi>w</mi><mi>m</mi></msub><mi>ψ</mi><mo stretchy="false">(</mo><msub><mi>r</mi><mi>i</mi></msub><mn>,</mn><msub><mi>E</mi><mi>g</mi></msub><mn>,</mn><msub><mi>Ω</mi><mi>m</mi></msub><mo stretchy="false">)</mo><mo>|</mo><msub><mi>λ</mi><mi>m</mi></msub><mo>|</mo><msub><mi>A</mi><mi>i</mi></msub></mrow><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>m</mi><mo>−</mo><mn>1</mn></mrow><mi>M</mi></munderover></mstyle><msub><mi>w</mi><mi>m</mi></msub><mi>ψ</mi><mo stretchy="false">(</mo><msub><mi>r</mi><mi>i</mi></msub><mn>,</mn><msub><mi>E</mi><mi>g</mi></msub><mn>,</mn><msub><mi>Ω</mi><mi>m</mi></msub><mo stretchy="false">)</mo><mo>|</mo><msub><mi>λ</mi><mi>m</mi></msub><mo>|</mo><msub><mi>A</mi><mi>i</mi></msub></mrow></mfrac><mn>,</mn></mrow></math> where ψ ( r i , E g , Ω m ) represents the angular flux in the spatial mesh i for energy group g and direction m . I is the total number of spatial mesh, G is the total number of energy groups, M is the total number of directions. λ m means direction cosine with respect to a set of orthogonal coordinate axes. w M is the corresponding directional weight. A i means the area of the spatial mesh i. The probability of any spatial mesh, p ( i ), is calculated by Eq. 1 based on the radio of the total number of particles crossing at that location to the total number of particles crossing the entire common surface. p ( g / i ) in Eq. 2 is the normalized energy spectrum. Eq. 3 defines the probability for every angle used in the discrete ordinates method. That is, if the calculation used an S8 quadrature scheme, then there would be 80 probabilities, each corresponding to a unique direction. <math display="block" id="M4"><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mrow><mi>j</mi><mo>−</mo><mn>1</mn></mrow></munderover></mstyle><msub><mi>p</mi><mi>i</mi></msub><mo>≤</mo><mi>ξ</mi><mo>≤</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>j</mi></munderover></mstyle><msub><mi>p</mi><mi>i</mi></msub><mn>.</mn></mrow></math> The space bin is chosen by Eq. 4, and ξ means a random number. The energy and angular bin selection were conducted in the same fashion as shown in Eq. 4, and the final energy and direction were sampled randomly from within the selected bin. III. PROGRAM SYSTEM 1 The flow chart of three-dimensional S N -MC coupled program system was shown in Fig. 1 Fig. 1 . The coupled program system consists of three-dimensional S N code, S N -MC interface program, MC source subroutine and MC code. In practical application, a common surface is first specified in both S N mesh model and the MC model. In S N calculation, the angular fluxes (as a function of mesh interval, energy group and discrete direction) on the common surface are stored in a BNDRYS formatted file. In MC simulation, this surface is specified as the source. The interface program reads the BNDRYS-format file generated by S N , and then uses the equation given in Section 2 to generate the necessary probability distributions. The MC source subroutine is linked to MC code to generate MC source particles variables by using the probability. The Cartesian ( X , Y , Z ) and Cylindrical ( R , θ , Z ) can be selected in the calculation by the coupling method. IV. CALCULATION RESULTS AND DISCUSSION 1 The HBR-2 benchmark was adopted to validate the capability of the 3D S N -MC coupling methodology. The specific activities for the six dosimeters were calculated. The dosimeters were irradiated during cycle 9 in the surveillance capsule. The surveillance capsule is located on the mid-place of the HBR-2 core. The experimental measurements and results of DORT code are included in the benchmark documentation. The 45 <math id="M5"><mrow><msup><mrow/><mo>∘</mo></msup></mrow></math> sector of 3D calculated model is displayed in Fig. 2 Fig. 2 . As shown in Fig. 2 Fig. 2 , the model was separated to S N mesh model and the MC model. The inner surface of core barrel is specified as the common surface. TORT [ 5 5 ] and MCNP [ 6 6 ] were used via the program system for S N and MC model calculation, respectively. The S N model was analyzed with the TEXT10 library, which have 30 neutron and 12 gamma energy groups. The point-wise cross-section files generated from ENDF/B-VI data library [ 7 7 ] were used in the MCNP calculations. For S N calculation, the order of scattering expansion is P3, and the angular quadrature set is S8. The mesh grids for the S N model are 110×36×48. The angular fluxes on the common surface, calculated by TORT, were stored in a BNDRY format file, and then the interface program and the source subroutine were employed to generate the source particles of MCNP. The MC model describes the detailed geometry and material compositions of surveillance capsule mounting on the thermal shield. The neutron flux values and activities at the dosimeter position of the surveillance capsule are calculated by S N -MC coupled method. It is seen from Fig. 3 Fig. 3 that the S N -MC neutron spectrum appears to be in relatively good agreement with the DORT-BUGLE96 solution in HBR-2 benchmark report. The ratios of the calculated specific activities by different codes with the measurements were shown in Table 1 Table 1 . The average ratios of calculated-to-measured ( C/M ) for the six dosimeters in the surveillance capsule for DORT, S N -MC and TORT are 0.90±0.04, 1.03±0.04 and 1.04±0.04, respectively. The C/M of DORT and TORT come from [ 4 4 ] and [ 8 8 ], respectively. Table 1 Table 1 shows that the S N -MC calculated specific activities are quite similar to the TORT results and the experimental measurements. As to be expected from the neutron flux comparisons, the S N -MC calculated specific activities are larger than the DORT results. The necessary simplification of certain features for the HBR-2 reactor, the different codes and cross-section library differences can all be the contributing factors. 1 1 Code 1 1 237 Np( n , f ) 137 Cs 1 1 238 U( n , f ) 137 Cs 1 1 58 Ni( n , p ) 58 Co 1 1 54 Fe( n , p ) 54 Mn 1 1 46 Ti( n , p ) 46 Sc 1 1 63 Cu( n , a ) 60 Co 1 1 Average 1 1 DORT 1 1 0.90 1 1 0.85 1 1 0.96 1 1 0.93 1 1 0.85 1 1 0.90 1 1 0.89±0.04 1 1 S N -MC 1 1 1.06 1 1 0.96 1 1 1.01 1 1 1.03 1 1 1.01 1 1 1.09 1 1 1.03±0.04 1 1 TORT 1 1 1.12 1 1 1.07 1 1 1.00 1 1 0.99 1 1 1.03 1 1 1.05 1 1 1.04±0.04 V. CONCLUSIONS 1 A program system has been developed to implement the coupled S N -MC scheme by integrating the 3D S N code TORT, the MC transport code MCNP, the coupling interface program and source subroutine. The program system is able to deal with the shielding problem with bulk shield and complex geometry in large dimensions. The HBR-2 benchmark problem has been performed to validate the program system. Good agreements were obtained by comparing with the measurements and the results of other codes. The successful application demonstrates that the coupling scheme with the program system is a powerful computational tool for three-dimensional shielding analysis of complex and large nuclear facilities. References Emmett M B , Burgart C E , Hoffman T J . ORNL-4853: DOMINO, a general purpose code for coupling discrete ordinates and Monte Carlo radiation transport calculations , USA , 1973 . ORNL-4853: DOMINO, a general purpose code for coupling discrete ordinates and Monte Carlo radiation transport calculations Eggleston J E , Abdou M A , Youssef M Z . Fusion Eng Des , 1998 , 42 : 281 - 288 . Fusion Eng Des Kurosawa M . Radiat Prot Dosim , 2005 , 116 : 513 - 517 . Radiat Prot Dosim Remec I , Kam F B K . ORNL/TM-13204: H.B. Robinson-2 pressure vessel benchmark , USA , 1997 . ORNL/TM-13204: H.B. Robinson-2 pressure vessel benchmark Rhoades W A , Simpson D B . ORNL/TM-13221: The TORT Three-dimensional Discrete Ordinates Neutron/Photon Transport Code , USA , 1997 . ORNL/TM-13221: The TORT Three-dimensional Discrete Ordinates Neutron/Photon Transport Code Briesmeister J F . LA-13709-M: MCNP:A general Monte Carlo N-particle transport code, version 4C , USA , 2000 . LA-13709-M: MCNP:A general Monte Carlo N-particle transport code, version 4C Wienke H , Herman M . IAEA-NDS-176: FENDL/MG-2.0 and FENDL/MC-2.0 , the processed cross-section libraries for neutron-photon transport calculations , version 1 , Austria , 1998 . IAEA-NDS-176: FENDL/MG-2.0 and FENDL/MC-2.0 Yang S H , Chen Y X , Wang W J , et al . Chin J Nucl Sci Eng , 2011 , 31 : 294 - 298 . Chin J Nucl Sci Eng 10.13538/j.1001-8042/nst.25.S010601 Supported by the CAP1400 Safety Evaluation Technology and Independent Verification Test of Ministry of Science and Technology Major Projects (2011ZX06002-010), CAP1400 Safety Review Key Technology Study of Ministry of Science and Technology Major Projects (2013ZX06002001). 2013-06-30 2014-09-15 2014-11-20 2014-12-25 2021-04-11 2014 Nuclear Science and Techniques S1 25