On-the-fly Doppler broadening method based on optimal double-exponential formula

On-the-fly Doppler broadening method based on optimal double-exponential formula 1 2 first-author 1 1 1 1 2 corresp Corresponding author Email: liqin.hu@fds.org.cn Corresponding author Email: liqin.hu@fds.org.cn liqin.hu@fds.org.cn Key Laboratory of Neutronics and Radiation Safety, Institute of Nuclear Energy Safety Technology, Chinese Academy of Sciences, Hefei 230031, China Key Laboratory of Neutronics and Radiation Safety, Institute of Nuclear Energy Safety Technology, Chinese Academy of Sciences, Hefei 230031, China University of Science and Technology of China, Hefei 230027, China University of Science and Technology of China, Hefei 230027, China As temperature changes constantly in nuclear reactor operation, on-the-fly Doppler broadening methods are commonly adopted for generating nuclear cross sections at various temperatures in neutron transport simulation. Among the existing methods, the widely-used SIGMA1 approach is inefficient because it involves error function and Taylor series expansion. In this paper, we present a new on-the-fly Doppler broadening with optimal double-exponential formula based on SuperMC to improve efficiency with given accuracy. In this method, double-exponential formula in 1/16 steps is used for broadening cross section at low energy, with both accuracy and efficiency. Meanwhile, the Gauss-Hermite quadrature of different orders is used for broadening cross section at resonance energy. The method can generate neutron cross section rapidly and precisely at the desired temperature. Typical nuclide cross sections and benchmarking tests are presented in detail. on-the-fly Doppler broadening cross section double-exponential formula SuperMC 1. Introduction 1 In nuclear and particle physics, neutron cross section is used to express the probability of interaction between the incident neutron and target nucleus, and it depends on the relative speed between them. Since the relative motion varies significantly due to thermal agitation of target nucleus, the resonance magnitude of cross section increases with the temperature, while the resonance peak magnitude decreases. This phenomenon, described as Doppler broadening, has a great influence on the total resonance absorption or fission in materials, causing intense effects of neutron transport as temperature changes. In operation of a nuclear reactor, temperature changes constantly in extensive ranges. For simulating neutronics behavior in a realistic reactor, considering the Doppler broadening effect, it is desirable to store all the related cross section data at any temperature, especially when thermal-hydraulic feedback is involved. A common approach is to pre-generate the cross section based on nuclear data processing systems such as NJOY[ 1 1 ] and PREPRO[ 2 2 ], which means tremendous memory footprint in the neutron transport simulation. To solve this problem, on-the-fly Doppler broadening methods have been developed to calculate cross section at any temperature in transport simulation, such as the SIGMA1 method [ 3 3 ], regression model method [ 4 4 ], target motion sampling method[ 5 5 ], multiple approximation method [ 6 6 ] and compact-extended model method [ 7 7 ]. The SIGMA1 is widely used due to its characteristic of low memory footprint. However, it is not efficient in evaluating complementary error function and Taylor series expansion. So, Gauss-Hermite quadrature is used to reduce the computation time of SIGMA1[ 8 8 ], but it is only suitable for broadening the cross section at resonance energy region. At low energy regions, the Gauss-Hermite quadrature will cause relative errors of up to 2% [ 9 9 ], which is much higher than the generally acceptable tolerance of 0.1%[ 10 10 ]. In this paper, a novel on-the-fly Doppler broadening method with optimal double-exponential formula is developed based on Super Monte Carlo Program for Nuclear and Radiation Simulation (SuperMC)[ 11 11 , 12 12 ], which is a general, intelligent, accurate and precise simulation software system for the nuclear design and safety evaluation of nuclear systems[ 13 13 - 15 15 13-15 ]. 2. Optimal double-exponential method 1 Generally, assuming the velocity of the target nuclei in the medium is isotropic Maxwellian distribution, the well-known Doppler broadening equation describes the effective cross section of the incident neutron <math id="M1"><mrow><msub><mi>σ</mi><mrow><msub><mi>T</mi><mn>2</mn></msub></mrow></msub><mrow><mo>(</mo><mi>v</mi><mo>)</mo></mrow></mrow></math> at speed <math id="M2"><mi>v</mi></math> and temperature T 2 from temperature T 1 as: <math display="block" id="M3"><mrow><msub><mi>σ</mi><mrow><msub><mi>T</mi><mn>2</mn></msub></mrow></msub><mrow><mo>(</mo><mi>v</mi><mo>)</mo></mrow><mo>=</mo><mo> </mo><mfrac><mn>1</mn><mrow><msup><mi>v</mi><mn>2</mn></msup></mrow></mfrac><msqrt><mrow><mfrac><mi>β</mi><mi>π</mi></mfrac></mrow></msqrt><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∫</mo></mstyle><mn>0</mn><mi>∞</mi></munderover><msubsup><mi>v</mi><mi>r</mi><mn>2</mn></msubsup><msub><mi>σ</mi><mrow><msub><mi>T</mi><mn>1</mn></msub></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>v</mi><mi>r</mi></msub></mrow><mo>)</mo></mrow><mrow><mo>[</mo><mrow><msup><mtext>e</mtext><mrow><mo>−</mo><mi>β</mi><msup><mrow><mrow><mo>(</mo><mrow><msub><mi>v</mi><mi>r</mi></msub><mo>−</mo><mi>v</mi></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow></msup><mo>−</mo><msup><mtext>e</mtext><mrow><mo>−</mo><mi>β</mi><msup><mrow><mrow><mo>(</mo><mrow><msub><mi>v</mi><mi>r</mi></msub><mo>+</mo><mi>v</mi></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow></msup></mrow><mo>]</mo></mrow><mtext>d</mtext><msub><mi>v</mi><mi>r</mi></msub><mo>,</mo></mrow></math> where v r is relative speed between incident neutron and target nuclide; and <math id="M4"><mrow><mtext>  </mtext><mi>β</mi><mo>=</mo><mi>A</mi><mo>/</mo><mrow><mo>[</mo><mrow><mi>k</mi><mrow><mo>(</mo><mrow><msub><mi>T</mi><mn>2</mn></msub><mo>−</mo><msub><mi>T</mi><mn>1</mn></msub></mrow><mo>)</mo></mrow></mrow><mo>]</mo></mrow></mrow></math> with <math id="M5"><mi>k</mi></math> being the Boltzmann’s constant, A=M/m being mass ratio between target ( M ) and neutron ( m ). <math id="M6"><mo> </mo></math> Defining <math id="M7"><mrow><msup><mi>x</mi><mn>2</mn></msup><mo>=</mo><mi>β</mi><msubsup><mi>v</mi><mi>r</mi><mn>2</mn></msubsup></mrow></math> and <math id="M8"><mrow><mtext> </mtext><msup><mi>y</mi><mn>2</mn></msup><mo>=</mo><mi>β</mi><msup><mi>v</mi><mn>2</mn></msup></mrow></math> , we have <math display="block" id="M9"><mrow><mtable><mtr><mtd><mrow><msub><mi>σ</mi><mrow><msub><mi>T</mi><mn>2</mn></msub></mrow></msub><mrow><mo>(</mo><mi>y</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mn>1</mn><mrow><msup><mi>y</mi><mn>2</mn></msup><msqrt><mi>π</mi></msqrt></mrow></mfrac><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∫</mo></mstyle><mn>0</mn><mi>∞</mi></munderover><msup><mi>x</mi><mn>2</mn></msup><msub><mi>σ</mi><mrow><msub><mi>T</mi><mn>1</mn></msub></mrow></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mrow><mo>[</mo><mrow><msup><mtext>e</mtext><mrow><mo>−</mo><msup><mrow><mrow><mo>(</mo><mrow><mi>x</mi><mo>−</mo><mi>y</mi></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow></msup><mo>−</mo><msup><mtext>e</mtext><mrow><mo>−</mo><msup><mrow><mrow><mo>(</mo><mrow><mi>x</mi><mo>+</mo><mi>y</mi></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow></msup></mrow><mo>]</mo></mrow><mtext>d</mtext><mi>x</mi></mrow></mtd></mtr></mtable><mo>.</mo></mrow></math> Eq.(2) can be simplified as <math id="M10"><mrow><msub><mi>σ</mi><mrow><msub><mi>T</mi><mn>2</mn></msub></mrow></msub><mrow><mo>(</mo><mi>y</mi><mo>)</mo></mrow><mo>=</mo><msubsup><mi>σ</mi><mrow><msub><mi>T</mi><mn>2</mn></msub></mrow><mo>*</mo></msubsup><mrow><mo>(</mo><mi>y</mi><mo>)</mo></mrow><mo>−</mo><msubsup><mi>σ</mi><mrow><msub><mi>T</mi><mn>2</mn></msub></mrow><mo>*</mo></msubsup><mrow><mo>(</mo><mrow><mo>−</mo><mi>y</mi></mrow><mo>)</mo></mrow></mrow></math> by defining <math display="block" id="M11"><mrow><msubsup><mi>σ</mi><mrow><msub><mi>T</mi><mn>2</mn></msub></mrow><mo>*</mo></msubsup><mrow><mo>(</mo><mi>y</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mn>1</mn><mrow><msup><mi>y</mi><mn>2</mn></msup></mrow></mfrac><msqrt><mrow><mfrac><mn>1</mn><mi>π</mi></mfrac></mrow></msqrt><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∫</mo></mstyle><mn>0</mn><mi>∞</mi></munderover><msup><mi>x</mi><mn>2</mn></msup><msub><mi>σ</mi><mrow><msub><mi>T</mi><mn>1</mn></msub></mrow></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><msup><mtext>e</mtext><mrow><mo>−</mo><msup><mrow><mrow><mo>(</mo><mrow><mi>x</mi><mo>−</mo><mi>y</mi></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow></msup><mtext>d</mtext><mi>x</mi><mo>.</mo></mrow></math> Cullen and Weisbin showed[ 16 16 ] that the range of integration is limited to y-4≤ x < y +4 for <math id="M12"><mrow><msubsup><mi>σ</mi><mrow><msub><mi>T</mi><mn>2</mn></msub></mrow><mo>*</mo></msubsup><mrow><mo>(</mo><mi>y</mi><mo>)</mo></mrow></mrow></math> and 0 < x < 4 for <math id="M13"><mrow><mtext> </mtext><msubsup><mi>σ</mi><mrow><msub><mi>T</mi><mn>2</mn></msub></mrow><mo>*</mo></msubsup><mrow><mo>(</mo><mrow><mo>−</mo><mi>y</mi></mrow><mo>)</mo></mrow></mrow></math> . For y > 4 (equivalent to E > 16 kT /A), <math id="M14"><mrow><msubsup><mi>σ</mi><mrow><msub><mi>T</mi><mn>2</mn></msub></mrow><mo>*</mo></msubsup><mrow><mo>(</mo><mrow><mo>−</mo><mi>y</mi></mrow><mo>)</mo></mrow></mrow></math> can be ignored due to integrand decays with the square of exponent. Above all, Eq. (2) can now be rewritten as <math display="block" id="M15"><mrow><msub><mi>σ</mi><mrow><msub><mi>T</mi><mn>2</mn></msub></mrow></msub><mrow><mo>(</mo><mi>y</mi><mo>)</mo></mrow><mo>=</mo><mrow><mo>{</mo><mrow><mtable><mtr><mtd><mrow><mfrac><mn>1</mn><mrow><msup><mi>y</mi><mn>2</mn></msup><msqrt><mi>π</mi></msqrt></mrow></mfrac><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∫</mo></mstyle><mn>0</mn><mrow><mi>y</mi><mo>+</mo><mn>4</mn></mrow></munderover><msup><mi>x</mi><mn>2</mn></msup><msub><mi>σ</mi><mrow><msub><mi>T</mi><mn>1</mn></msub></mrow></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mrow><mo>[</mo><mrow><msup><mtext>e</mtext><mrow><mo>−</mo><msup><mrow><mrow><mo>(</mo><mrow><mi>x</mi><mo>−</mo><mi>y</mi></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow></msup><mo>−</mo><msup><mtext>e</mtext><mrow><mo>−</mo><msup><mrow><mrow><mo>(</mo><mrow><mi>x</mi><mo>+</mo><mi>y</mi></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow></msup></mrow><mo>]</mo></mrow><mtext>d</mtext><mi>x</mi><mo>,</mo><mi>y</mi><mo>≤</mo><mn>4</mn></mrow></mtd></mtr><mtr><mtd><mrow><mfrac><mn>1</mn><mrow><msup><mi>y</mi><mn>2</mn></msup><msqrt><mi>π</mi></msqrt></mrow></mfrac><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∫</mo></mstyle><mn>0</mn><mi>∞</mi></munderover><msup><mi>x</mi><mn>2</mn></msup><msub><mi>σ</mi><mrow><msub><mi>T</mi><mn>1</mn></msub></mrow></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><msup><mtext>e</mtext><mrow><mo>−</mo><msup><mrow><mrow><mo>(</mo><mrow><mi>x</mi><mo>−</mo><mi>y</mi></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow></msup><mtext>d</mtext><mi>x</mi><mo>,</mo><mi>y</mi><mo>></mo><mn>4</mn></mrow></mtd></mtr></mtable></mrow></mrow><mo>.</mo></mrow></math> The SIGMA1 method is inefficient because the analytical form of Eq. (4) introduces complementary error function and Taylor series expansion to avoid the loss of accuracy. In order to improve the efficiency of SIGMA1, a novel on-the-fly Doppler broadening method with optimal double-exponential (DE) formula is proposed for efficiency problem. In the method, concerning the large error occurs at low energy (y≤4) of Gauss-Hermite quadrature, double-exponential formula with optimal step size was proposed for both accuracy and efficiency. Meanwhile, Gauss-Hermite quadrature with different orders was optimized for broadening cross section at resonance energy (y>4). 2.1 Low energy treatment 2 As mentioned above, the DE formula is proposed for broadening low energy cross section, which is used to solve the integrals with high accuracy [ 17 17 - 18 18 17-18 ]. For the integration <math id="M16"><mrow><mtext> </mtext><mi>I</mi><mo>=</mo><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∫</mo></mstyle><mi>a</mi><mi>b</mi></munderover><mi>f</mi><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mtext>d</mtext><mi>t</mi></mrow></math> , the formula applies a variable substitution with <math id="M17"><mrow><mi>t</mi><mo>=</mo><mi>ψ</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow></mrow></math> as: <math display="block" id="M18"><mrow><mi>I</mi><mo>=</mo><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∫</mo></mstyle><mi>a</mi><mi>b</mi></munderover><mi>f</mi><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mtext>d</mtext><mi>t</mi><mo>=</mo><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∫</mo></mstyle><mrow><mo>−</mo><mi>∞</mi></mrow><mi>∞</mi></munderover><mi>f</mi><mrow><mo>(</mo><mrow><mi>ψ</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow></mrow><mo>)</mo></mrow><mi>ψ</mi><mo>'</mo><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mtext>d</mtext><mi>x</mi><mo>,</mo></mrow></math> where the transformation <math id="M19"><mrow><mi>ψ</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow></mrow></math> and <math id="M20"><mrow><mi>ψ</mi><mo>'</mo><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow></mrow></math> of DE formula are defined by <math display="block" id="M21"><mrow><mtable><mtr><mtd><mrow><mi>ψ</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mrow><mi>b</mi><mo>−</mo><mi>a</mi></mrow><mn>2</mn></mfrac><mi>tanh</mi><mrow><mo>(</mo><mrow><mfrac><mtext>π</mtext><mn>2</mn></mfrac><mi>sinh</mi><mi>x</mi></mrow><mo>)</mo></mrow><mo>+</mo><mfrac><mrow><mi>b</mi><mo>+</mo><mi>a</mi></mrow><mn>2</mn></mfrac></mrow></mtd></mtr><mtr><mtd><mrow><mi>ψ</mi><mo>'</mo><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mrow><mrow><mo>(</mo><mrow><mi>b</mi><mo>−</mo><mi>a</mi></mrow><mo>)</mo></mrow></mrow><mn>2</mn></mfrac><mfrac><mrow><mtext>π</mtext><mi>cosh</mi><mi>x</mi></mrow><mrow><mn>2</mn><msup><mrow><mi>cosh</mi></mrow><mn>2</mn></msup><mrow><mo>(</mo><mrow><mfrac><mi>π</mi><mn>2</mn></mfrac><mi>sinh</mi><mi>x</mi></mrow><mo>)</mo></mrow></mrow></mfrac></mrow></mtd></mtr></mtable><mo>.</mo></mrow></math> The trapezoidal formula with N and an equal step size h are applied to solve Eq. (5) as: <math display="block" id="M22"><mrow><mtable><mtr><mtd><mrow><msub><mi>I</mi><mi>h</mi></msub><mo>=</mo><mi>h</mi><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>k</mi><mo>=</mo><mo>−</mo><mi>N</mi><mo>/</mo><mn>2</mn></mrow><mrow><mi>N</mi><mo>/</mo><mn>2</mn></mrow></munderover><mi>f</mi><mrow><mo>(</mo><mrow><mi>ψ</mi><mrow><mo>(</mo><mrow><mi>k</mi><mi>h</mi></mrow><mo>)</mo></mrow></mrow><mo>)</mo></mrow><mi>ψ</mi><mo>'</mo><mrow><mo>(</mo><mrow><mi>k</mi><mi>h</mi></mrow><mo>)</mo></mrow></mrow></mtd></mtr></mtable><mo>.</mo></mrow></math> For ease of use in implementation , we rewrite DE formula by defining the abscissas x k and the weights w k as: <math display="block" id="M23"><mrow><mtable><mtr><mtd><mrow><msub><mi>x</mi><mi>k</mi></msub><mo>=</mo><mtext>tanh</mtext><mo stretchy="false">(</mo><mfrac><mtext>π</mtext><mn>2</mn></mfrac><mi>sinh</mi><mi>k</mi><mi>h</mi><mo stretchy="false">)</mo></mrow></mtd></mtr><mtr><mtd><mrow><msub><mi>w</mi><mi>k</mi></msub><mo>=</mo><mfrac><mrow><mfrac><mtext>π</mtext><mn>2</mn></mfrac><mi>cosh</mi><mi>k</mi><mi>h</mi></mrow><mrow><msup><mrow><mi>cosh</mi></mrow><mn>2</mn></msup><mrow><mo>(</mo><mrow><mfrac><mtext>π</mtext><mn>2</mn></mfrac><mi>sinh</mi><mi>k</mi><mi>h</mi></mrow><mo>)</mo></mrow></mrow></mfrac></mrow></mtd></mtr></mtable></mrow></math> Accordingly, Eq. (7) can be rewritten as <math display="block" id="M24"><mrow><msub><mi>I</mi><mi>h</mi></msub><mo>=</mo><mfrac><mrow><mi>b</mi><mo>−</mo><mi>a</mi></mrow><mn>2</mn></mfrac><mi>h</mi><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>k</mi><mo>=</mo><mo>−</mo><mi>N</mi><mo>/</mo><mn>2</mn></mrow><mrow><mi>N</mi><mo>/</mo><mn>2</mn></mrow></munderover><mi>f</mi><mrow><mo>(</mo><mrow><mfrac><mrow><mi>b</mi><mo>−</mo><mi>a</mi></mrow><mn>2</mn></mfrac><msub><mi>x</mi><mi>k</mi></msub><mo>+</mo><mfrac><mrow><mi>b</mi><mo>+</mo><mi>a</mi></mrow><mn>2</mn></mfrac></mrow><mo>)</mo></mrow><msub><mi>w</mi><mi>k</mi></msub><mo>.</mo></mrow></math> As shown in Fig. 1 Fig. 1 , the definitional domain of x k and w k can be shrank to [-3,3] for the most requirement due to rapid convergence and accuracy. For on-the-fly Doppler broadening purpose, according to the formula and range of the integration, Eq. (4) at <math id="M25"><mrow><mtext>y</mtext><mo>≤</mo><mn>4</mn><mtext> </mtext></mrow></math> can be expressed as <math display="block" id="M26"><mrow><mtable><mtr><mtd><mrow><mo> </mo><mo> </mo><msub><mi>σ</mi><mrow><msub><mi>T</mi><mn>2</mn></msub><mo>,</mo><mi>y</mi><mo>≤</mo><mn>4</mn></mrow></msub><mrow><mo>(</mo><mi>y</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mn>1</mn><mrow><msup><mi>y</mi><mn>2</mn></msup><msqrt><mi>π</mi></msqrt></mrow></mfrac><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∫</mo></mstyle><mn>0</mn><mrow><mi>y</mi><mo>+</mo><mn>4</mn></mrow></munderover><msup><mi>x</mi><mn>2</mn></msup><msub><mi>σ</mi><mrow><msub><mi>T</mi><mn>1</mn></msub></mrow></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mrow><mo>[</mo><mrow><msup><mtext>e</mtext><mrow><mo>−</mo><msup><mrow><mrow><mo>(</mo><mrow><mi>x</mi><mo>−</mo><mi>y</mi></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow></msup><mo>−</mo><msup><mtext>e</mtext><mrow><mo>−</mo><msup><mrow><mrow><mo>(</mo><mrow><mi>x</mi><mo>+</mo><mi>y</mi></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow></msup></mrow><mo>]</mo></mrow><mtext>d</mtext><mi>x</mi></mrow></mtd></mtr><mtr><mtd><mrow><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo>≈</mo><mfrac><mn>2</mn><mrow><msup><mi>y</mi><mn>2</mn></msup><msqrt><mi>π</mi></msqrt></mrow></mfrac><mi>h</mi><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>k</mi><mo>=</mo><mo>−</mo><mi>N</mi><mo>/</mo><mn>2</mn></mrow><mrow><mi>N</mi><mo>/</mo><mn>2</mn></mrow></munderover><mi>f</mi><mrow><mo>(</mo><mrow><mfrac><mrow><mi>y</mi><mo>+</mo><mn>4</mn></mrow><mn>2</mn></mfrac><msub><mi>x</mi><mi>k</mi></msub><mo>+</mo><mfrac><mrow><mi>y</mi><mo>+</mo><mn>4</mn></mrow><mn>2</mn></mfrac></mrow><mo>)</mo></mrow><msub><mi>w</mi><mi>k</mi></msub><mo> </mo></mrow></mtd></mtr></mtable><mo>,</mo></mrow></math> where <math id="M27"><mrow><mtext> </mtext><mi>f</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>=</mo><msup><mi>x</mi><mn>2</mn></msup><msub><mi>σ</mi><mrow><msub><mi>T</mi><mn>1</mn></msub></mrow></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mrow><mo>[</mo><mrow><msup><mtext>e</mtext><mrow><mo>−</mo><msup><mrow><mrow><mo>(</mo><mrow><mi>x</mi><mo>−</mo><mi>y</mi></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow></msup><mo>−</mo><msup><mtext>e</mtext><mrow><mo>−</mo><msup><mrow><mrow><mo>(</mo><mrow><mi>x</mi><mo>+</mo><mi>y</mi></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow></msup></mrow><mo>]</mo></mrow></mrow></math> . As the tabulated low energy cross section ends before reaching the lower boundary <math id="M28"><mrow><mo stretchy="false">[</mo><mrow><mo>(</mo><mrow><mi>y</mi><mo>+</mo><mn>4</mn></mrow><mo>)</mo></mrow><mo>/</mo><mn>2</mn><msub><mi>x</mi><mrow><mtext>min</mtext></mrow></msub><mo>+</mo><mrow><mo>(</mo><mrow><mi>y</mi><mo>+</mo><mn>4</mn></mrow><mo>)</mo></mrow><mo>/</mo><mn>2</mn><mo stretchy="false">]</mo></mrow></math> , the cross section is approximated as 1/ v up to zero, and Eq. (10) becomes <math display="block" id="M29"><mrow><mtable><mtr><mtd><mrow><msub><mi>σ</mi><mrow><msub><mi>T</mi><mn>2</mn></msub><mo>,</mo><mi>y</mi><mo>≤</mo><mn>4</mn></mrow></msub><mrow><mo>(</mo><mi>y</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mn>2</mn><mrow><msup><mi>y</mi><mn>2</mn></msup><msqrt><mtext>π</mtext></msqrt></mrow></mfrac><mi>h</mi><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>k</mi><mo>=</mo><mo>−</mo><mi>N</mi><mo>/</mo><mn>2</mn></mrow><mrow><mi>N</mi><mo>/</mo><mn>2</mn></mrow></munderover><mi>f</mi><mrow><mo>(</mo><mrow><mfrac><mrow><mi>y</mi><mo>+</mo><mn>4</mn></mrow><mn>2</mn></mfrac><msub><mi>x</mi><mrow><mtext>min</mtext></mrow></msub><mo>+</mo><mfrac><mrow><mi>y</mi><mo>+</mo><mn>4</mn></mrow><mn>2</mn></mfrac></mrow><mo>)</mo></mrow><msub><mi>w</mi><mi>k</mi></msub></mrow></mtd></mtr></mtable></mrow></math> In order to analysis influences arose by changes of the step sizes to the efficiency and accuracy of DE method, Fig. 2 Fig. 2 shows the maximum relative error and broadened efficiency ratio between DE and SIGMA1 method in generating total cross section of 1 H, 208 Pb and 238 U with diminishing step sizes of h= 1/8, 1/16, 1/32, 1/64 and 1/128 at temperatures of 600, 900 and 1200 K. The broadened efficiency ratio here is defined as <math id="M30"><mrow><mtext>broadened efficiency ratio</mtext><mo>=</mo><mfrac><mrow><mtext>broadening time of SIMGA</mtext><mn>1</mn><mtext> </mtext><mrow><mo>(</mo><mrow><mtext>y</mtext><mo>≤</mo><mn>4</mn></mrow><mo>)</mo></mrow></mrow><mrow><mtext>broadening time of DE </mtext><mrow><mo>(</mo><mrow><mtext>y</mtext><mo>≤</mo><mn>4</mn></mrow><mo>)</mo></mrow></mrow></mfrac><mo>.</mo></mrow></math> From Fig. 2 Fig. 2 , as the interval increases, the accuracy improves but the broadened efficiency ratio decline rapidly. Within the acceptable tolerance value of 0.1%, step size of 1/16 is chosen as the optimal value for both efficiency and accuracy. To be more specific, Table 1 Table 1 lists the simulated maximum relative errors for broadening total cross section with h =1/16 for 1 H, 10 B, 16 O, 90 Zr, 208 Pb, 232 Th, 235 U and 238 U, at 600, 900 and 1200 K, which proves the accuracy of the method. 1 1 Nuclides 1 1 600 K 1 1 900 K 1 1 1200 K 1 1 1 H 1 1 0.00264 1 1 0.00128 1 1 0.00178 1 1 10 B 1 1 0.00582 1 1 0.00475 1 1 0.00359 1 1 16 O 1 1 0.00173 1 1 0.00176 1 1 0.00243 1 1 90 Zr 1 1 0.00245 1 1 0.00155 1 1 0.0013 1 1 208 Pb 1 1 0.00434 1 1 0.00187 1 1 0.00154 1 1 232 Th 1 1 0.00126 1 1 0.0038 1 1 0.00217 1 1 235 U 1 1 0.00355 1 1 0.00576 1 1 0.00382 1 1 238 U 1 1 0.00251 1 1 0.00263 1 1 0.00248 2.2 Resonance energy treatment 2 The Gauss-Hermite quadrature is used for broadening cross section at resonance energy region. Due to the limitation of this method, it can be used for broadening cross section of y>4 only. In this case, Eq. (4) can be rewrite as: <math display="block" id="M31"><mrow><msub><mi>σ</mi><mrow><msub><mi>T</mi><mn>2</mn></msub><mo>,</mo><mi>y</mi><mo>></mo><mn>4</mn></mrow></msub><mrow><mo>(</mo><mi>y</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mn>1</mn><mrow><msup><mi>y</mi><mn>2</mn></msup><msqrt><mi>π</mi></msqrt></mrow></mfrac><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>k</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><msub><mi>w</mi><mi>k</mi></msub><msup><mrow><mrow><mo>(</mo><mrow><msub><mi>z</mi><mi>k</mi></msub><mo>+</mo><mi>y</mi></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup><msub><mi>σ</mi><mrow><msub><mi>T</mi><mn>1</mn></msub></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>z</mi><mi>k</mi></msub><mo>+</mo><mi>y</mi></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> Where Z k is the abscissas of the Hermite polynomial H n ( z ) with n degrees and <math id="M32"><mrow><msub><mi>w</mi><mi>k</mi></msub><mo>=</mo><msup><mn>2</mn><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msup><mi>n</mi><mo>!</mo><msqrt><mi>π</mi></msqrt><mo>/</mo><mrow><mo>[</mo><mrow><msup><mi>n</mi><mn>2</mn></msup><msub><mi>H</mi><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>z</mi><mi>k</mi></msub></mrow><mo>)</mo></mrow></mrow><mo>]</mo></mrow></mrow></math> . Since the orders selection of Gauss-Hermite quadrature affects greatly the efficiency and accuracy, Fig. 3 Fig. 3 shows the maximum relative error and broadened efficiency ratio in generating total cross section of 16 O and 238 U with different Gauss-Hermite orders at 600, 900 and 1200 K. In this section, the broadened efficiency ratio is defined as <math id="M33"><mrow><mtext>broadened efficiency ratio</mtext><mo>=</mo><mfrac><mrow><mtext>broadening time of SIMGA</mtext><mn>1</mn><mtext> </mtext><mo stretchy="false">(</mo><mi>y</mi><mo>></mo><mn>4</mn><mo stretchy="false">)</mo></mrow><mrow><mtext>broadening time of Gauss</mtext><mo>−</mo><mtext>Hermit </mtext><mo stretchy="false">(</mo><mtext>y</mtext><mo>></mo><mn>4</mn><mo stretchy="false">)</mo></mrow></mfrac><mo>.</mo></mrow></math> In Fig.3 Fig.3 , as the orders increases, the accuracy improves but the broadening efficiency declines rapidly. Within the acceptable tolerance value of 0.1%, different orders of Gauss-Hermite quadrature can be used in the premise of efficiency. The order of 6 is chosen for 16 O with optimal performance, whereas the order of 16 is optimal for 238 U for better efficiency. Based on the above analysis, orders of Gauss-Hermite quadrature are analyzed and determined for each nuclide. Table 2 Table 2 gives maximum relative error for broadening total cross section with different optimal orders of typical nuclides at 600, 900 and 1200 K. Comparing to constant order selecting, it can greatly improve the efficiency for broadening cross section of different nuclides with optimal orders. 1 1 Nuclides 1 1 Orders 1 1 600 K 1 1 900 K 1 1 1200 K 1 1 1 H 1 1 2 1 1 0.00120 1 1 0.00037 1 1 0.00076 1 1 10 B 1 1 2 1 1 0.03833 1 1 0.03856 1 1 0.03555 1 1 11 Na 1 1 4 1 1 0.01799 1 1 0.04032 1 1 0.09207 1 1 16 O 1 1 6 1 1 0.03988 1 1 0.05624 1 1 0.06884 1 1 35 Cl 1 1 14 1 1 0.08001 1 1 0.06128 1 1 0.09800 1 1 70 Ge 1 1 4 1 1 0.00121 1 1 0.00448 1 1 0.00230 1 1 90 Zr 1 1 20 1 1 0.08411 1 1 0.09214 1 1 0.09878 1 1 232 Th 1 1 15 1 1 0.04059 1 1 0.05019 1 1 0.09140 1 1 234 U 1 1 16 1 1 0.04303 1 1 0.05077 1 1 0.09014 1 1 235 U 1 1 15 1 1 0.03658 1 1 0.04187 1 1 0.08502 1 1 238 U 1 1 16 1 1 0.04671 1 1 0.05436 1 1 0.08591 1 1 239 Pu 1 1 16 1 1 0.03674 1 1 0.04546 1 1 0.07699 3. Results and Discussion 1 The accuracy and efficiency of the DE method are verified with the cross section broadening of typical nuclides, criticality safety benchmarks and Doppler coefficient of reactivity, in terms of nuclear data and calculation results of transport. 3.1 Doppler broadening of typical nuclides cross sections tests 2 The nuclides of 10 B, 16 O, 56 Fe, 90 Zr, 232 Th, 234 U, 235 U, 238 U and 239 Pu are used to compare the accuracy and efficiency of Doppler broadening cross sections between DE and SIGMA1 method. Fig.4 Fig.4 shows the broadened relative error of 238 U at 300–600, 900 and 1200 K. Although it increases gradually with the neutron energy due to computational precision, the relative error is always less than the acceptable tolerance of 0.1%. The broadened efficiency ratio between the DE and SIGMA1 method is shown in Fig. 5 Fig. 5 . The low energy region and resonance region are used for efficiency comparison, and the broadened efficiency ratio is defined as <math id="M34"><mrow><mtext>Broadened efficiency ratio</mtext><mo>=</mo><mfrac><mrow><mtext>Broadening time of SIMGA</mtext><mn>1</mn><mtext> </mtext><mrow><mo>(</mo><mrow><mi>y</mi><mo>∈</mo><mrow><mo>(</mo><mrow><mn>0</mn><mo>,</mo><mo>+</mo><mi>∞</mi></mrow><mo>)</mo></mrow></mrow><mo>)</mo></mrow></mrow><mrow><mtext>Broadening time of DE </mtext><mrow><mo>(</mo><mrow><mtext>y</mtext><mo>≤</mo><mn>4</mn></mrow><mo>)</mo></mrow><mtext> and Gauss</mtext><mo>−</mo><mtext>Hermit </mtext><mo stretchy="false">(</mo><mtext>y</mtext><mo>></mo><mn>4</mn><mo stretchy="false">)</mo></mrow></mfrac><mo>.</mo></mrow></math> In Fig.5 Fig.5 , average efficiency of the DE method was improved by a factor of 12. For time consumption comparison such as 238 U at 1200 and 1800 K on an Intel Core i5–4590 processor, the broadening times of total cross section with SIGMA1 are 3.965 and 5.192 s, whereas those with the DE method are only 0.269 and 0.275 s, respectively. 3.2 Criticality safety benchmarks 2 For evaluating the on-the-fly Doppler broadening method, eight criticality safety benchmark cases from ICSBEP[ 19 19 ] are selected to testify the accuracy. From 300 K to 600 K, the cross sections of each nuclide are broadened, and each case has 500 cycles with 100,000 histories per cycle. In the calculation, the first 50 inactive cycles are skipped due to statistical error. Table 3 Table 3 shows the k eff results of the DE and SIGMA1 methods. All results are in the statistical margin of deviations and the maximum deviation is 15 pcm. It proves the validity and reliability of the DE method. 1 1 Cases 1 1   1 1 DE 1 1 SIGMA1 1 1 Dev. 1 1 PU-MET-FAST 1 1 001 1 1 0.99974±0.00009 1 1 0.99968±0.00008 1 1 6 pcm 1 1   1 1 006 1 1 1.00123±0.00011 1 1 1.00122±0.00011 1 1 1 pcm 1 1   1 1 009 1 1 1.00517±0.00010 1 1 1.00506±0.00010 1 1 11 pcm 1 1 HEU-MET-FAST 1 1 001 1 1 1.00001±0.00009 1 1 0.99995±0.00009 1 1 6 pcm 1 1   1 1 003 1 1 0.99548±0.00009 1 1 0.99546±0.00009 1 1 9 pcm 1 1   1 1 008 1 1 0.99602±0.00009 1 1 0.99591±0.00009 1 1 11 pcm 1 1 U233-MET-FAST 1 1 001 1 1 0.99954±0.00009 1 1 0.99957±0.00009 1 1 3 pcm 1 1   1 1 002 1 1 0.99920±0.00009 1 1 0.99905±0.00009 1 1 15 pcm 3.3 Doppler coefficient of reactivity 2 The benchmark of Doppler coefficient of reactivity[ 20 20 ] is shown in Fig. 6 Fig. 6 . The temperature of fuel, cladding and moderator is 600 K at hot zero power (HZP) condition. However, the fuel temperature is 900 K at hot full power (HFP) condition. The Doppler coefficient is simply defined as <math display="block" id="M35"><mrow><mi>D</mi><mi>C</mi><mo>=</mo><mfrac><mrow><mi>Δ</mi><msub><mi>ρ</mi><mrow><mtext>Dop</mtext></mrow></msub></mrow><mrow><mi>Δ</mi><msub><mi>T</mi><mrow><mtext>fuel</mtext></mrow></msub></mrow></mfrac><mo>,</mo></mrow></math> where <math id="M36"><mrow><mo> </mo><mi>Δ</mi><msub><mi>ρ</mi><mrow><mtext>Dop</mtext></mrow></msub><mo>=</mo><mrow><mo>(</mo><mrow><msub><mi>k</mi><mrow><mtext>HFP</mtext></mrow></msub><mo>−</mo><msub><mi>k</mi><mrow><mtext>HZP</mtext></mrow></msub></mrow><mo>)</mo></mrow><mo>/</mo><mrow><mo>(</mo><mrow><msub><mi>k</mi><mrow><mtext>HFP</mtext></mrow></msub><mtext> </mtext><msub><mi>k</mi><mrow><mtext>HZP</mtext></mrow></msub></mrow><mo>)</mo></mrow></mrow></math> and <math id="M37"><mrow><mi>Δ</mi><msub><mi>T</mi><mrow><mtext>fuel</mtext></mrow></msub><mo>=</mo></mrow></math> 300 K. According to uncertain factors of the benchmark in safety analyses, the uncertainty of measurement is less than 10%. Each of the cases has 500 cycles with 20,000 histories per cycle. In the calculation, the first 50 inactive cycles are skipped due to statistical error. Fig. 7 Fig. 7 shows the comparison for Doppler coefficient of reactivity with DE method and reference value. The Doppler coefficient results agree with the reference values in maximum relative error of < 3.74%, well meeting the accuracy requirements of the benchmark model. 4. Conclusion 1 We have proposed a novel on-the-fly Doppler broadening method with optimal double-exponential formula based on SuperMC, which uses the double-exponential formula in 1/16 steps at low energy and the Gauss-Hermite quadrature with different orders at resonance energy. For validation the method, cross section of typical nuclides, criticality safety benchmarks from ICSBEP as well as Doppler coefficient of reactivity are used. The results indicate that the method can generate neutron cross section rapidly and precisely at desired temperature. Within the acceptable tolerance value of 0.1%, the average broadened efficiency of the method is improved by a factor of 12. The results prove the method can satisfy the need of on-the-fly Doppler broadening. References R. Macfarlane , D. Muir , R. Boicourt et al ., The NJOY nuclear data processing system, version 2012, LA-UR-12-27079 , Los Alamos National Laboratory ( 2012 ). The NJOY nuclear data processing system, version 2012, LA-UR-12-27079 D. Cullen , PREPRO 2015: 2015 ENDF/B Pre-processing Codes, IAEA-NDS-39 , Rev. 14 International Atomic Energy Agency ( 2015 ). PREPRO 2015: 2015 ENDF/B Pre-processing Codes, IAEA-NDS-39 D. Cullen , Program SIGMA1 (version 79-1): Doppler broaden evaluated cross sections in the evaluated nuclear data file/version B (ENDF/B) format, UCRL-50400 , Lawrence Livermore National Laboratory ( 1979 ). Program SIGMA1 (version 79-1): Doppler broaden evaluated cross sections in the evaluated nuclear data file/version B (ENDF/B) format, UCRL-50400 G. Yesilyurt , W. Martin , F. Brown , On-the-fly Doppler broadening for Monte Carlo codes . Nucl. Sci. Eng . 171 , 239 - 257 ( 2012 ). doi: 10.13182/NSE11-67 http://doi.org/10.13182/NSE11-67 On-the-fly Doppler broadening for Monte Carlo codes T. Viitanen , J. Leppanen , Explicit Treatment of thermal motion in continuous-energy Monte Carlo tracking routines . Nucl. Sci. Eng . 171 , 165 - 173 ( 2012 ). doi: 10.13182/NSE11-36 http://doi.org/10.13182/NSE11-36 Explicit Treatment of thermal motion in continuous-energy Monte Carlo tracking routines B. Forget , S. Xu , K. Smith , Direct Doppler broadening in Monte Carlo simulations using the multipole representation . Ann. Nucl. Energy 64 , 78 - 85 ( 2014 ). doi: 10.1016/j.anucene.2013.09.043 http://doi.org/10.1016/j.anucene.2013.09.043 Direct Doppler broadening in Monte Carlo simulations using the multipole representation A. Villanueva , J. Granada , Experimental evidence sustaining a compact-extended model for Doppler broadening of neutron absorption resonances . Ann. Nucl. Energy 38 , 1389 - 1398 ( 2011 ). doi: 10.1016/j.anucene.2011.01.034 http://doi.org/10.1016/j.anucene.2011.01.034 Experimental evidence sustaining a compact-extended model for Doppler broadening of neutron absorption resonances P. Romano , T. Trumbull , Comparison of algorithms for Doppler broadening pointwise tabulated cross section . Ann. Nucl. Energy 75 , 358 - 364 ( 2015 ). doi: 10.1016/j.anucene.2014.08.046 http://doi.org/10.1016/j.anucene.2014.08.046 Comparison of algorithms for Doppler broadening pointwise tabulated cross section C. Dean , R. Perry , R. Neal et al ., Validation of run-time Doppler broadening in MONK with JEFF3.1 . J. Korean Phys. Soc . 59 , 1162 - 1165 ( 2011 ). doi: 10.3938/jkps.59.1162 http://doi.org/10.3938/jkps.59.1162 Validation of run-time Doppler broadening in MONK with JEFF3.1 G. Ferran , W. Haeck , A new method for the Doppler Broadening of the Solbrig’s kernel using a Fourier transforms . Nucl. Sci. Eng . 179 , 1 - 17 ( 2015 ). doi: 10.13182/NSE14-64 http://doi.org/10.13182/NSE14-64 A new method for the Doppler Broadening of the Solbrig’s kernel using a Fourier transforms Y. Wu , J. Song , H. Zheng et al ., CAD-based Monte Carlo program for integrated simulation of nuclear system SuperMC . Ann. Nucl. Energy 82 , 161 - 168 ( 2015 ). doi: 10.1016/j.anucene.2014.08.058 http://doi.org/10.1016/j.anucene.2014.08.058 CAD-based Monte Carlo program for integrated simulation of nuclear system SuperMC Y. Wu , FDS Team , CAD-based interface programs for fusion neutron transport simulation . Fusion Eng. Des . 84 , 1987 - 1992 ( 2009 ). doi: 10.1016/j.fusengdes.2008.12.041 http://doi.org/10.1016/j.fusengdes.2008.12.041 CAD-based interface programs for fusion neutron transport simulation Y. Wu , Z. Chen , L. Hu et al ., Identification of safety gaps for fusion demonstration reactors . Nature Energy 1 , 16154 ( 2016 ). doi: 10.1038/nenergy.2016.154 http://doi.org/10.1038/nenergy.2016.154 Identification of safety gaps for fusion demonstration reactors Y. Wu , J. Jiang , M. Wang et al ., A fusion-driven subcritical system concept based on viable technologies . Nucl. Fusion 51 , 103036 ( 2011 ). doi: 10.1088/0029-5515/51/10/103036 http://doi.org/10.1088/0029-5515/51/10/103036 A fusion-driven subcritical system concept based on viable technologies Y. Wu , FDS Team , Design Analysis of the China Dual-Functional Lithium Lead (DFLL) test blanket module in ITER . Fusion Eng. Des . 82 , 1893 - 1903 ( 2007 ). doi: 10.1016/j.fusengdes.2007.08.012 http://doi.org/10.1016/j.fusengdes.2007.08.012 Design Analysis of the China Dual-Functional Lithium Lead (DFLL) test blanket module in ITER D. Cullen , C. Weisbin , Exact Doppler broadening of tabulated cross section . Nucl. Sci. Eng . 60 , 199 - 229 ( 1976 ). doi: 10.13182/NSE76-1 http://doi.org/10.13182/NSE76-1 . Exact Doppler broadening of tabulated cross section H. Takahasi , M. Mori , Double-exponential formulas for numerical integration . Publ. Res. I. Math. Sci . 9 , 721 - 741 ( 1974 ). Double-exponential formulas for numerical integration M. Mori , M. Sugihara , The double-exponential transformation in numerical analysis . J. Comput. Appl. Math . 127 , 287 - 296 ( 2001 ). doi: 10.1016/S0377-0427(00)00501-X http://doi.org/10.1016/S0377-0427(00)00501-X The double-exponential transformation in numerical analysis J. Briggs , L. Scott , A. Nouri , The International Criticality Safety Benchmark Evaluation Project . Nucl. Sci. Eng . 145 , 1 - 10 ( 2003 ). doi: 10.13182/NSE03-14 http://doi.org/10.13182/NSE03-14 The International Criticality Safety Benchmark Evaluation Project R. Mosteller , ENDF/B-V, ENDF/B-VI, ENDF/B-VII.0 Results for the Doppler-Defect Benchmark, LA-UR-07-0922 , Los Alamos National Laboratory ( 2007 ). ENDF/B-V, ENDF/B-VI, ENDF/B-VII.0 Results for the Doppler-Defect Benchmark, LA-UR-07-0922 The authors would like to show their great appreciation to other members of the FDS Team for supports to this research. 10.1007/s41365-017-0318-4 This work was supported by the Strategic Priority Science & Technology Program of the Chinese Academy of Sciences (No. XDA03040000), the Innovation Foundation of the Chinese Academy of Sciences (No. CXJJ-16Q231), the National Natural Science Foundation of China (NSFC) (Nos. 11305205, 11305203, 11405204 and 11605233), the National Magnetic Confinement Fusion Science Program of China (No. 2014GB112001), the Special Program for Informatization of the Chinese Academy of Sciences (No. XXH12504-1-09), the Anhui Provincial Special project for High Technology Industry, and the Special Project of Youth Innovation Promotion Association of Chinese Academy of Sciences, the Industrialization Fund. 2017-01-10 2017-02-22 2017-03-02 2017-10-30 2017-11-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-0318-4 http://dx.doi.org/10.1007/s41365-017-0318-4 . 2017 Nuclear Science and Techniques 11 28