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

NUCLEAR ENERGY SCIENCE AND ENGINEERING

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

Rui Chen，

Li-Juan Hao，

Bin Wu，

Jing Song，

Li-Qin Hu

Nuclear Science and Techniques

Vol.28, No.11

Article number 166

Published in print 01 Nov 2017

Available online 30 Oct 2017

DOI：10.1007/s41365-017-0318-4

37500

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 broadeningcross sectiondouble-exponential formulaSuperMC

1. Introduction

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] and PREPRO[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], regression model method [4], target motion sampling method[5], multiple approximation method [6] and compact-extended model method [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], 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], which is much higher than the generally acceptable tolerance of 0.1%[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,12], which is a general, intelligent, accurate and precise simulation software system for the nuclear design and safety evaluation of nuclear systems[13-15].

2. Optimal double-exponential method

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 $σ_{T_{2}} (v)$ at speed $v$ and temperature T₂ from temperature T₁ as:

σ_{T_{2}} (v) = \frac{1}{v^{2}} \sqrt{\frac{β}{π}} \int_{0}^{\infty} v_{r}^{2} σ_{T_{1}} (v_{r}) [e^{- β {(v_{r} - v)}^{2}} - e^{- β {(v_{r} + v)}^{2}}] d v_{r},

(1)

where v_r is relative speed between incident neutron and target nuclide; and $β = A / [k (T_{2} - T_{1})]$ with $k$ being the Boltzmann’s constant, A=M/m being mass ratio between target (M) and neutron (m). Defining $x^{2} = β v_{r}^{2}$ and $y^{2} = β v^{2}$ , we have

\begin{matrix} σ_{T_{2}} (y) = \frac{1}{y^{2} \sqrt{π}} \int_{0}^{\infty} x^{2} σ_{T_{1}} (x) [e^{- {(x - y)}^{2}} - e^{- {(x + y)}^{2}}] d x \end{matrix} .

(2)

Eq.(2) can be simplified as $σ_{T_{2}} (y) = σ_{T_{2}}^{*} (y) - σ_{T_{2}}^{*} (- y)$ by defining

σ_{T_{2}}^{*} (y) = \frac{1}{y^{2}} \sqrt{\frac{1}{π}} \int_{0}^{\infty} x^{2} σ_{T_{1}} (x) e^{- {(x - y)}^{2}} d x .

(3)

Cullen and Weisbin showed[16] that the range of integration is limited to y-4≤x<y+4 for $σ_{T_{2}}^{*} (y)$ and 0<x<4 for $σ_{T_{2}}^{*} (- y)$ . For y>4 (equivalent to E>16kT/A), $σ_{T_{2}}^{*} (- y)$ can be ignored due to integrand decays with the square of exponent. Above all, Eq. (2) can now be rewritten as

σ_{T_{2}} (y) = {\begin{matrix} \frac{1}{y^{2} \sqrt{π}} \int_{0}^{y + 4} x^{2} σ_{T_{1}} (x) [e^{- {(x - y)}^{2}} - e^{- {(x + y)}^{2}}] d x, y \leq 4 \\ \frac{1}{y^{2} \sqrt{π}} \int_{0}^{\infty} x^{2} σ_{T_{1}} (x) e^{- {(x - y)}^{2}} d x, y > 4 \end{matrix} .

(4)

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

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-18].

For the integration $I = \int_{a}^{b} f (t) d t$ , the formula applies a variable substitution with $t = ψ (x)$ as:

I = \int_{a}^{b} f (t) d t = \int_{- \infty}^{\infty} f (ψ (x)) ψ' (x) d x,

(5)

where the transformation $ψ (x)$ and $ψ' (x)$ of DE formula are defined by

\begin{matrix} ψ (x) = \frac{b - a}{2} \tanh (\frac{π}{2} \sinh x) + \frac{b + a}{2} \\ ψ' (x) = \frac{(b - a)}{2} \frac{π \cosh x}{2 \cosh^{2} (\frac{π}{2} \sinh x)} \end{matrix} .

(6)

The trapezoidal formula with N and an equal step size h are applied to solve Eq. (5) as:

\begin{matrix} I_{h} = h \sum_{k = - N / 2}^{N / 2} f (ψ (k h)) ψ' (k h) \end{matrix} .

(7)

For ease of use in implementation, we rewrite DE formula by defining the abscissas x_k and the weights w_k as:

\begin{matrix} x_{k} = tanh (\frac{π}{2} \sinh k h) \\ w_{k} = \frac{\frac{π}{2} \cosh k h}{\cosh^{2} (\frac{π}{2} \sinh k h)} \end{matrix}

(8)

Accordingly, Eq. (7) can be rewritten as

I_{h} = \frac{b - a}{2} h \sum_{k = - N / 2}^{N / 2} f (\frac{b - a}{2} x_{k} + \frac{b + a}{2}) w_{k} .

(9)

As shown in 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.

Fig. 1

The variation tendency of x_k and w_k

For on-the-fly Doppler broadening purpose, according to the formula and range of the integration, Eq. (4) at $y \leq 4$ can be expressed as

\begin{matrix} σ_{T_{2}, y \leq 4} (y) = \frac{1}{y^{2} \sqrt{π}} \int_{0}^{y + 4} x^{2} σ_{T_{1}} (x) [e^{- {(x - y)}^{2}} - e^{- {(x + y)}^{2}}] d x \\ \approx \frac{2}{y^{2} \sqrt{π}} h \sum_{k = - N / 2}^{N / 2} f (\frac{y + 4}{2} x_{k} + \frac{y + 4}{2}) w_{k} \end{matrix},

(10)

where $f (x) = x^{2} σ_{T_{1}} (x) [e^{- {(x - y)}^{2}} - e^{- {(x + y)}^{2}}]$ . As the tabulated low energy cross section ends before reaching the lower boundary $[(y + 4) / 2 x_{min} + (y + 4) / 2]$ , the cross section is approximated as 1/v up to zero, and Eq. (10) becomes

\begin{matrix} σ_{T_{2}, y \leq 4} (y) = \frac{2}{y^{2} \sqrt{π}} h \sum_{k = - N / 2}^{N / 2} f (\frac{y + 4}{2} x_{min} + \frac{y + 4}{2}) w_{k} \end{matrix}

(11)

In order to analysis influences arose by changes of the step sizes to the efficiency and accuracy of DE method, Fig. 2 shows the maximum relative error and broadened efficiency ratio between DE and SIGMA1 method in generating total cross section of ¹H, ²⁰⁸Pb and ²³⁸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

Fig. 2.

Maximum relative error and broadened efficiency ratio in generating total cross section of ¹H (a), ²⁰⁸Pb (b) and ²³⁸U (c) with diminishing intervals (h=1/8, 1/16, 1/32,1/64 and 1/128) at 600, 900 and 1200 K.

$broadened efficiency ratio = \frac{broadening time of SIMGA 1 (y \leq 4)}{broadening time of DE (y \leq 4)} .$

From 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 lists the simulated maximum relative errors for broadening total cross section with h=1/16 for ¹H, ¹⁰B, ¹⁶O, ⁹⁰Zr, ²⁰⁸Pb, ²³²Th, ²³⁵U and ²³⁸U, at 600, 900 and 1200 K, which proves the accuracy of the method.

Maximum relative error for broadened total cross section of ¹H, ¹⁰B, ¹⁶O, ⁹⁰Zr, ²⁰⁸Pb, ²³²Th, ²³⁵U and ²³⁸U, at h=1/16 of the DE formula.

Nuclides	600 K	900 K	1200 K
¹H	0.00264	0.00128	0.00178
¹⁰B	0.00582	0.00475	0.00359
¹⁶O	0.00173	0.00176	0.00243
⁹⁰Zr	0.00245	0.00155	0.0013
²⁰⁸Pb	0.00434	0.00187	0.00154
²³²Th	0.00126	0.0038	0.00217
²³⁵U	0.00355	0.00576	0.00382
²³⁸U	0.00251	0.00263	0.00248

2.2 Resonance energy treatment

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:

σ_{T_{2}, y > 4} (y) = \frac{1}{y^{2} \sqrt{π}} \sum_{k = 1}^{N} w_{k} {(z_{k} + y)}^{2} σ_{T_{1}} (z_{k} + y),

(12)

Where Z_k is the abscissas of the Hermite polynomial H_n(z) with n degrees and $w_{k} = 2^{n - 1} n! \sqrt{π} / [n^{2} H_{n - 1} (z_{k})]$ . Since the orders selection of Gauss-Hermite quadrature affects greatly the efficiency and accuracy, Fig. 3 shows the maximum relative error and broadened efficiency ratio in generating total cross section of ¹⁶O and ²³⁸U with different Gauss-Hermite orders at 600, 900 and 1200 K. In this section, the broadened efficiency ratio is defined as

Fig. 3.

Maxrelative error (solid line) and broadened efficiency ratio (dashed line) in generating total cross section of ¹⁶O (a) and ²³⁸U (b) with different Gauss-Hermite orders at 600 K, 900 K and 1200 K.

$broadened efficiency ratio = \frac{broadening time of SIMGA 1 (y > 4)}{broadening time of Gauss - Hermit (y > 4)} .$

In 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 ¹⁶O with optimal performance, whereas the order of 16 is optimal for ²³⁸U for better efficiency.

Based on the above analysis, orders of Gauss-Hermite quadrature are analyzed and determined for each nuclide. 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.

Maximum relative error for broadened cross section of the typical nuclides at different orders of Gauss quadrature.

Nuclides	Orders	600 K	900 K	1200 K
¹H	2	0.00120	0.00037	0.00076
¹⁰B	2	0.03833	0.03856	0.03555
¹¹Na	4	0.01799	0.04032	0.09207
¹⁶O	6	0.03988	0.05624	0.06884
³⁵Cl	14	0.08001	0.06128	0.09800
⁷⁰Ge	4	0.00121	0.00448	0.00230
⁹⁰Zr	20	0.08411	0.09214	0.09878
²³²Th	15	0.04059	0.05019	0.09140
²³⁴U	16	0.04303	0.05077	0.09014
²³⁵U	15	0.03658	0.04187	0.08502
²³⁸U	16	0.04671	0.05436	0.08591
²³⁹Pu	16	0.03674	0.04546	0.07699

3. Results and Discussion

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

The nuclides of ¹⁰B, ¹⁶O, ⁵⁶Fe, ⁹⁰Zr, ²³²Th, ²³⁴U, ²³⁵U, ²³⁸U and ²³⁹Pu are used to compare the accuracy and efficiency of Doppler broadening cross sections between DE and SIGMA1 method.

Fig.4 shows the broadened relative error of ²³⁸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%.

Fig. 4.

Relativeerrors of broadened total cross section of ²³⁸U at 300–600 K (a), 900 K (b) and 1200 K (c), by DE and SIGMA1 methods.

The broadened efficiency ratio between the DE and SIGMA1 method is shown in Fig. 5. The low energy region and resonance region are used for efficiency comparison, and the broadened efficiency ratio is defined as

Fig. 5

Thebroadened efficiency ratio between DE and SIGMA1 method at 300–600, 900, 1200, 1800 and 3000 K.

$Broadened efficiency ratio = \frac{Broadening time of SIMGA 1 (y \in (0, + \infty))}{Broadening time of DE (y \leq 4) and Gauss - Hermit (y > 4)} .$

In Fig.5, average efficiency of the DE method was improved by a factor of 12. For time consumption comparison such as ²³⁸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

For evaluating the on-the-fly Doppler broadening method, eight criticality safety benchmark cases from ICSBEP[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 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.

The k_eff results and deviations of ICSBEP cases with DE and SIGMA1 methods

Cases		DE	SIGMA1	Dev.
PU-MET-FAST	001	0.99974±0.00009	0.99968±0.00008	6 pcm
	006	1.00123±0.00011	1.00122±0.00011	1 pcm
	009	1.00517±0.00010	1.00506±0.00010	11 pcm
HEU-MET-FAST	001	1.00001±0.00009	0.99995±0.00009	6 pcm
	003	0.99548±0.00009	0.99546±0.00009	9 pcm
	008	0.99602±0.00009	0.99591±0.00009	11 pcm
U233-MET-FAST	001	0.99954±0.00009	0.99957±0.00009	3 pcm
	002	0.99920±0.00009	0.99905±0.00009	15 pcm

3.3 Doppler coefficient of reactivity

The benchmark of Doppler coefficient of reactivity[20] is shown in 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

Fig. 6

Model of Doppler coefficient of reactivity.

D C = \frac{Δ ρ_{Dop}}{Δ T_{fuel}},

(12)

where $Δ ρ_{Dop} = (k_{HFP} - k_{HZP}) / (k_{HFP} k_{HZP})$ and $Δ T_{fuel} =$ 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 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.

Fig. 7

Doppler coefficient for normal and enriched UO₂ fuel

4. Conclusion

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).

D. Cullen, PREPRO 2015: 2015 ENDF/B Pre-processing Codes, IAEA-NDS-39, Rev. 14 International Atomic Energy Agency (2015).

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).

G. Yesilyurt, W. Martin, F. Brown,

On-the-ﬂy Doppler broadening for Monte Carlo codes

. Nucl. Sci. Eng. 171, 239-257(2012). doi: 10.13182/NSE11-67

Baidu Scholar

Google Scholar

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

Baidu Scholar

Google Scholar

B. Forget, S. Xu, K. Smith,

Direct Doppler broadening in Monte Carlo simulations using the multipole representation