SP3-coupled global variance reduction method based on RMC code

NUCLEAR ENERGY SCIENCE AND ENGINEERING

SP3-coupled global variance reduction method based on RMC code

Qing-Quan Pan ，

Teng-Fei Zhang，

Xiao-Jing Liu，

Hui He，

Kan Wang

Nuclear Science and Techniques

Vol.32, No.11

Article number 122

Published in print 01 Nov 2021

Available online 16 Nov 2021

DOI：10.1007/s41365-021-00973-0

35900

A global variance reduction (GVR) method based on the SPN method is proposed. First, global multi-group cross-sections are obtained by Monte Carlo (MC) global homogenization. Then, the SP3 equation is solved to obtain the global flux distribution. Finally, the global weight windows are approximated by the global flux distribution, and the GVR simulation is performed. This GVR method is implemented as an automatic process in the RMC code. The SP3-coupled GVR method was tested on a modified version of the C5G7 benchmark with a thickened water shield. The results show that the SP3-coupled GVR method can improve the efficiency of the MC criticality calculation.

Video Abstract

RMC codeGlobal homogenizationVariance reductionSPN theory

INTRODUCTION

The reactor criticality calculation is an essential step in the physical design of a reactor. The reactor burnup and neutronic/thermal (N-TH) coupling calculations [1] are based on the reactor criticality calculation. If the accuracy of the criticality calculation cannot be guaranteed, the accuracy of the reactor burnup and N-TH coupling calculations will also not be guaranteed. The Monte Carlo (MC) method is widely used in reactor criticality calculations because of its high fidelity. However, a non-uniform distribution of the relative errors has been found in MC criticality calculations. This non-uniform distribution of the relative errors has a negative influence on the efficiency and accuracy of some local tallies. In some cases, it even causes numerical instability in the MC multi-physics simulation [2] and MC kinetic calculation [3]. To improve the efficiency and accuracy of MC criticality calculations and ensure the numerical stability of MC multi-physics simulations, global variance reduction (GVR) methods should be studied.

There are many variance reduction (VR) methods, such as Zheng’s method [4, 5], the adaptive method [6], the FW-CADIS method [7], the automated weight windows method [8], and the importance function method [9]. These methods are well-known for their value in VR applications. However, these methods were initially developed for MC fixed-source calculations and not for MC criticality calculations. Some VR methods have also been developed for MC criticality calculations. They include the UFS [10], UTD [11], and UVM [12] methods. These methods are easily implemented and have good universality. However, the UFS method assumes that the fission neutrons are uniformly distributed in the regions with fuel and have little influence on the non-fuel regions. Moreover, only the source bias technique and not the weight window technique is adopted in the UFS method. This limits its value for GVR. The UTD method can reduce the variance of various statistics, but because its basic principle is still the same as that of the UFS method, it faces the same limitations. The UVM method is newly proposed, and can improve the efficiency of the MC criticality calculation and flatten the global variance. However, the development of this method is still in the initial stages and needs to be further studied.

A SP3-coupled GVR method is proposed here to improve the accuracy and efficiency of the MC criticality calculation. The SP3-coupled method can be used for GVR in the MC criticality calculation and allows the analysis of the reactor from its interior to its exterior through a single MC criticality calculation. Using the SP3-coupled GVR method, an accurate radiation source can be provided for reactor shielding design through the MC criticality. Moreover, the SP3-coupled GVR method can be coupled with the reactor burnup calculation to analyze the influence of reactor dynamics on reactor shielding.

The SP3 method is a numerical method for solving the neutron transport equation with a simplified treatment of spatial anisotropy, and has been widely studied [13, 14]. The traditional SP3 method starts from a one-dimensional model. The one-dimensional governing equation is first found and the one-dimensional operator is subsequently replaced with a three-dimensional operator to obtain the three-dimensional governing equation. Compared with traditional transport theory, the spatial anisotropy is simplified in the SP3 method to achieve higher calculation efficiency. Compared with traditional diffusion theory, the SP3 method provides a better treatment of spatial anisotropy to achieve higher calculation accuracy. Therefore, the SP3 method balances accuracy and efficiency, and has great potential for coupling with the MC method to achieve accelerated calculations and GVR.

In this paper, a hybrid MC method based on SP3 theory is studied and a GVR method driven by the SP3 theory is developed for the RMC code [15]. Unlike the traditional two-step reactor physical calculation method, a one-step process of global MC homogenization and global reactor calculation is achieved in this method, and an automatic calculation framework for GVR is constructed. Unlike other hybrid MC methods, such as the fission matrix [16] and response matrix [17] methods, the neutron transport equation is solved instead of a matrix relationship. This provides advantages in terms of numerical stability, development potential, and application scope. The remainder of this paper is organized as follows: Section 2 introduces the theory, and provides the related formulas and details. Section 3 presents the numerical results with some basic analysis. A discussion is given in Section 4 before the conclusion.

METHODS

The SP3-coupled GVR method is implemented as an automatic process in the RMC code. This implies that only the input card of the RMC code is required, and the GVR calculation will subsequently be automatically completed by the RMC code. The flowchart for the calculation is shown in Figure 1. The calculation consists of three main steps:

Fig. 1

Flowchart of the SP3-coupled GVR method.

(1) The global cross-sections are determined by MC global homogenization.

(2) The SP3 equation is solved to obtain the global flux distribution.

(3) The global weight windows are determined by the global flux.

Each of the three steps is described in the following subsections.

2.1

MC Global Homogenization

The SP3-coupled GVR method is a hybrid MC method in which the MC method is coupled with the SP3 method. A deterministic reactor code [18] for solving the SP3 equation, NLSP3, is developed. In contrast to the discrete ordinate method (SN), global multi-group cross-sections (MGCS) are required to solve the SP3 equation. Therefore, the MGCS must first be determined. MC global homogenization was adopted to determine the MGCS for the following reasons:

(1) Because a continuous-energy nuclear library is used, no resonance effect treatment [19] is required.

(2) Because the homogenization process is not based on an assembly model, no leakage correction [20] is required.

(3) Because the geometric model is not approximated, a higher accuracy can be achieved compared with traditional assembly homogenization.

(4) Because the global homogenization is based on the entire reactor, the cross-sections of all the material regions can be determined through a single calculation.

In particular, Point (4) implies that MC global homogenization is suitable for an automated process. Because the MGCS used for the MC global homogenization have a significant influence on the global flux distribution, a high accuracy of the MGCS must be ensured. Meanwhile, because the cross-sections in all the material regions are determined simultaneously, a non-uniform distribution of the relative errors will lead to different cross-section precisions in different material regions. To improve the accuracy of the cross-sections and ensure that the cross-sections in all material regions have the same precision, the following measures are taken [21]:

(1)A fixed-source calculation with a flat neutron source is performed for MC global homogenization. The spatial distribution of neutrons is positively correlated with the neutron flux distribution if the VR technique is adopted in the MC simulation. In a large reactor, the neutron fluxes in different phase spaces may differ, especially when the reflective and shielding layers are considered. These differences result in a variation in the spatial distribution of neutrons. Therefore, although performing global homogenization is a reasonable approach in the MC criticality calculation, it is inevitable that the cross-sections in some regions will have higher precision compared to other regions.

The spatial distribution of neutrons can be artificially set in a fixed-source calculation. Suppose a flat neutron source is provided in the fixed-source calculation. In this case, all the material regions have the same density of neutron histories, and the cross sections obtained will all have the same precision. Therefore, the MC global homogenization in the fixed-source calculation is performed with a flat neutron source.

(2)The neutron fission reaction is restricted during neutron transport. Various reactions such as scattering, fission, and absorption occur during neutron transport. The reactor becomes critical when the loss of a neutron because of the absorption of the neutron or its transportation outside the system is compensated by the corresponding production of a neutron in a fission reaction. Typically, a fission reaction produces two or three neutrons.

In a supercritical system, the loss of one neutron is accompanied by the production of more than one neutron. Therefore, if fission neutrons are produced during the fixed-source calculation, the number of simulated neutrons will increase, and the calculation will never be terminated. As stated earlier, global homogenization is performed during the MC fixed-source calculation. To ensure the feasibility of the fixed-source calculation, neutron fission is restricted during MC global homogenization: The neutrons only perform random walks in the system, and no fission reaction occurs. However, the various reaction probabilities are recorded during neutron transport.

(3)The critical spectrum is approximated using the Watt fission spectrum. The MC homogenization involves a process of compressing and grouping the neutron spectra according to the energy. The neutron spectrum is essential in MC global homogenization. However, because the critical spectrum is unknown before the fixed-source calculation, it can only be approximated.

As stated before, neutron fission is restricted during the neutron transport process. This restriction leads to errors. Therefore, the neutron spectrum can be artificially adjusted according to the fission spectrum when the flat neutron source is set during the fixed-source calculation. It is known that the fission spectrum of U-235 is described by the Watt spectrum [22]. Therefore, the external fission source is distributed according to the Watt spectrum

p (E) = C e^{- E / a} \sinh \sqrt{b E},

(1)

where C, a, and b are coefficients that are ordinarily equal to 0.453 MeV, 0.965 MeV, and 2.29 MeV^-1, respectively.

It should be noted that although the Watt spectrum may not be the best one for the external source and is not universal for different fission nuclides, it is the best approximation that we can come up with.

(4) Corrections are adopted to make the neutron spectrum asymptotic to the criticality spectrum. The main consideration in the Watt fission spectrum are the fission reactions that occur during the neutron transport process. The Watt fission spectrum is a good approximation in the reactor core region where many fission reactions occur. However, in regions without fissile materials, the neutron spectrum is mainly influenced by scattering reactions, and the error will be larger when the Watt spectrum is used to approximate the critical spectrum. Therefore, there should be as many scattering reactions as possible before the neutrons are lost. After many collisions, the neutron energy spectrum in the reflective and shielding layers become close to the critical spectrum. To make the neutron spectrum asymptotically approach the critical spectrum, the cutoff weight limit is reduced when more scattering reactions occur. In the SP3-coupled GVR method in this study, the weight cutoff threshold was set to 1.0E-8, which is the default value in the RMC code for the VR simulation. In the nonanalog MC method, the scattering reactions are effectively increased by reducing the cutoff weight windows. After multiple collisions in the reflective and shielding layers, the energy spectrum of the neutrons becomes close to the critical spectrum. However, at this stage, the weight of the neutrons is small. The neutrons therefore have little influence on the tallied results. To increase the influence of the neutrons after multiple collisions on the cross sections, the reaction rate is recorded based on the number of reactions rather than the weight of the neutrons.

We calculated the MGCS using a fixed-source calculation, which is an approximation. Thus, the MGCS obtained inherently has errors, which implies that the efficiency is improved at the cost of compromised calculation accuracy. The SP3 solution only provides the parameters for VR, and the final result is given by the MC criticality calculation. Thus, even if the accuracy of the MGCS solution is compromised, the error in the MGCS solution does not affect the accuracy of the final result. This is one of the reasons why we adopt an approximate method to solve for the MGCS.

The SP3 equation is solved to obtain the global flux distribution. The total grouped cross-sections $Σ_{t, g}$ , grouped absorption cross-sections $Σ_{a, g}$ , grouped fission production cross-sections $ν Σ_{f, g}$ , inter-group scattering matrix $Σ_{s, g^{'} \to g}$ , and fission spectrum $χ_{g}$ need to be input to the NLSP3 code. Therefore, only these five grouped cross-sections are calculated in the MC global homogenization. It should be stated that it is not easy to obtain the anisotropic cross-sections during the MC global homogenization. To reduce the computational cost, only isotropic cross-sections are considered.

The grouped flux distribution and grouped reaction rates are tallied to calculate the grouped cross-sections. The traditional volume-weighted method is employed in which the MGCS is defined as

{\bar{Σ}}_{x, g}^{i} \equiv \frac{\int_{V_{i}} Σ_{x, g} ϕ_{g} (r) d V}{\int_{V_{i}} ϕ_{g} (r) d V} x = a, f, s, \cdot \cdot \cdot g = 1, \cdot \cdot \cdot, G,

(2)

where r is a spatial variable, φ_g(r) is the global flux distribution, and the subscript i in V_i labels the different spatial regions.

The track length is used to tally the flux distribution and reaction rates. The flux is calculated as

ϕ_{g} = \frac{\int_{E_{g}}^{E_{g - 1}} d E \int_{V} d V \sum_{i = 1}^{N} W \cdot T L_{V}^{i} (E)}{V \sum_{i = 1}^{N} W_{0}^{i}},

(3)

and the reaction rate as

R_{g} = Σ_{g} ϕ_{g} = \frac{\int_{E_{g}}^{E_{g - 1}} d E \int_{V} d V \sum_{i = 1}^{N} W \cdot T L_{V}^{i} (E) \cdot Σ (r, E)}{V \sum_{i = 1}^{N} W_{0}^{i}} .

(4)

The homogenized cross-section is therefore

Σ_{g} = \frac{\int_{E_{g}}^{E_{g - 1}} d E \int_{V} d V \sum_{i = 1}^{N} W \cdot T L_{V}^{i} (E) \cdot Σ (r, E)}{\int_{E_{g}}^{E_{g - 1}} d E \int_{V} d V \sum_{i = 1}^{N} W \cdot T L_{V}^{i} (E)},

(5)

where E_g and E_g-1 are the lower and upper bounds of the g^th energy interval, respectively; $W \cdot T L_{V}^{i}$ is the product of the weight and track length of the i^th neutron in region V; and N is the total number of simulated neutrons.

The calculation of the fission energy spectrum is relatively simple. The fission energy spectrum is calculated with reference to the number of fission neutrons produced in each energy group. The number of fission neutrons produced is then normalized to obtain the energy-dependent probability of neutron production in the energy interval E. The fission energy spectrum is equal to the energy-dependent probability of neutron production, that is,

χ (E) = p (E),

(6)

where p(E) is the energy-dependent probability of neutron production.

The process of calculating the fission matrix is the same as that for calculating the fission energy spectrum. The neutron energies before and after scattering are tallied to calculate the scattering probability:

P_{g \to g^{'}} = \frac{\int_{E_{g^{'}}}^{E_{g^{'} - 1}} d E^{'} \int_{E_{g}}^{E_{g - 1}} d E \int_{V} ϕ (r, E) Σ_{s} (r, E \to E^{'}) d V}{\int_{E_{g}}^{E_{g - 1}} d E \int_{V} ϕ (r, E) Σ_{s} (r, E) d V} .

(7)

Using the scattering probability, the inter-group transfer cross-section is calculated as

Σ_{s, g \to g^{'}} = P_{g \to g^{'}} \cdot Σ_{s, g} .

(8)

The MGCS is calculated using Equations (3)–(8) during the MC global homogenization. The MGCS of all material regions required by the NLSP3 code can thus be determined simultaneously in a single fixed-source calculation for subsequent use in the SP3 global calculation.

2.2

Solution of the SP3 Equation

The SP3 equation is solved to obtain the global flux distribution. The SP3 equation is adopted because of the following reasons:

(1) Compared with the diffusion and CMFD equations, the SP3 equation has higher accuracy [23].

(2) Solving the SP3 equation takes approximately the same amount of time as solving the diffusion equation [24].

(3) The SP3 equation is easier to apply to larger systems because it requires less memory and computing time compared to the SN equation.

(4) To date, there has been no research on coupling the SP3 theory and the MC method for GVR. For academic interest, we would like to make an attempt to investigate this topic.

The SPN theory [13, 14] is a simplified spherical harmonics method for solving the neutron transport equation. It is widely studied because it has great potential for next-generation reactor physics calculations. There are two variants of the SPN theory, namely, the traditional and the rigorous SPN theories. In the traditional SPN theory, the one-dimensional governing equation is first found from a one-dimensional model, and the one-dimensional operator is subsequently replaced with a three-dimensional operator to obtain the three-dimensional governing equation. The traditional SPN theory lacks a rigorous mathematical basis. The rigorous SPN theory is derived based on the variational method and the concept of grouped angular flux. The angular flux is clearly defined in the rigorous SPN theory so that the physical fundamentals are completely considered in the treatment of the boundary conditions. These two SPN theories share the same SPN equation but use different boundary conditions.

Although a complete mathematical derivation is provided in the rigorous SPN, the solution of the rigorous SPN equation suffers from numerical instability because of its incompatibility with the traditional transverse integral technique [25]. The traditional SPN equation is therefore solved in the SP3-coupled GVR method. The “N” in “SPN” represents the order of the polynomial expansion. When N is equal to 3, the SP3 equations are obtained as:

- D_{0, g}^{k} \nabla^{2} ϕ_{0, g}^{k} (r) + Σ_{r 0, g}^{k} ϕ_{0, g}^{k} (r) - 2 Σ_{r 0, g}^{k} ϕ_{2, g}^{k} (r) = S_{0, g}^{k} (r),

(9a)

- D_{2, g}^{k} \nabla^{2} ϕ_{2, g}^{k} (r) + Σ_{r 2, g}^{k} ϕ_{0, g}^{k} (r) - \frac{2}{5} Σ_{r 0, g}^{k} ϕ_{0, g}^{k} (r) = - \frac{2}{5} S_{0, g}^{k} (r),

(9b)

where the superscript k is the nodal index, the subscript g is the energy-interval index, $ϕ_{0, g}^{k}$ is the zeroth-order scalar flux at the k^th node in the g^th energy interval, $ϕ_{2, g}^{k}$ is the second-order scalar flux at the k^th node in the g^th energy interval, and

D_{0, g}^{k} = \frac{1}{3 Σ_{t, g}^{k}},

(10a)

D_{2, g}^{k} = \frac{9}{35 Σ_{t, g}^{k}},

(10b)

Σ_{r 0, g}^{k} = Σ_{t, g}^{k} - Σ_{s, g g^{'}}^{k},

(10c)

Σ_{r 2, g}^{k} = \frac{9}{5} Σ_{t, g}^{k} - \frac{4}{5} Σ_{s, g g^{'}}^{k} .

(10d)

The source is

S_{0, g}^{k} (r) = Σ_{\begin{array}{l} g^{'} = 1 \\ g^{'} \neq g \end{array}}^{G} Σ_{0, g^{'} g}^{k} [ϕ_{0, g}^{k} (r) - 2 ϕ_{2, g}^{k} (r)] + \frac{1}{K_{ef f}} χ_{g}^{k} Σ_{g^{'} = 1}^{G} ν Σ_{f, g^{'}}^{k} [ϕ_{0, g^{'}}^{k} (r) - 2 ϕ_{2, g^{'}}^{k} (r)] .

(11)

The boundary conditions [26] for the traditional SP3 equation are

J_{0}^{\pm} (r) = \frac{1}{4} ϕ_{0}^{M} (r) \pm \frac{1}{2} \overset{⇀}{n} \cdot J_{0}^{M} (r) - \frac{3}{16} ϕ_{2}^{M} (r),

(12a)

J_{2}^{\pm} (r) = \frac{21}{80} ϕ_{2}^{M} (r) \pm \frac{1}{2} \overset{⇀}{n} \cdot J_{2}^{M} (r) - \frac{3}{80} ϕ_{0}^{M} (r),

(12b)

where the superscript M represents the quantities at the node boundary, and $J_{0}^{\pm}$ and $J_{2}^{\pm}$ are the zeroth-order and second-order partial currents, respectively.

The numerical instability problem was encountered when the SP3 equation was solved using the traditional nonlinear iterative method [27]. The zeroth-order flux is usually hundreds of times larger than the second-order flux. When the traditional coupling corrective relationship was used to calculate the corrective factors, there was also a hundred-fold difference between the zeroth-order and second-order corrective factors, which led to an ill-conditioned iteration matrix and convergence failure. The following new coupling corrective relationship is proposed to solve the numerical instability problem :

\begin{array}{l} J_{0, g u +}^{k} = - D_{1, g u +}^{k, FDM} (f_{0, g u -}^{k + 1} {\bar{ϕ}}_{0, g}^{k + 1} - f_{0, g u +}^{k} {\bar{ϕ}}_{0, g}^{k}) - D_{1, g u +}^{k, N O D} (f_{0, g u -}^{k + 1} {\bar{ϕ}}_{0, g}^{k + 1} + f_{0, g u +}^{k} {\bar{ϕ}}_{0, g}^{k}) \\ - D_{3, g u +}^{k, FDM} (f_{2, g u -}^{k + 1} {\bar{ϕ}}_{2, g}^{k + 1} - f_{2, g u +}^{k} {\bar{ϕ}}_{2, g}^{k}) - D_{3, g u +}^{k, N O D} (f_{2, g u -}^{k + 1} {\bar{ϕ}}_{2, g}^{k + 1} + f_{2, g u +}^{k} {\bar{ϕ}}_{2, g}^{k}), \end{array}

(13a)

\begin{array}{l} J_{2, g u +}^{k} = - D_{2, g u +}^{k, FDM} (f_{2, g u -}^{k + 1} {\bar{ϕ}}_{2, g}^{k + 1} - f_{2, g u +}^{k} {\bar{ϕ}}_{2, g}^{k}) - D_{2, g u +}^{k, N O D} (f_{2, g u -}^{k + 1} {\bar{ϕ}}_{2, g}^{k + 1} + f_{2, g u +}^{k} {\bar{ϕ}}_{2, g}^{k}) \\ - D_{4, g u +}^{k, FDM} (f_{0, g u -}^{k + 1} {\bar{ϕ}}_{0, g}^{k + 1} - f_{0, g u +}^{k} {\bar{ϕ}}_{0, g}^{k}) - D_{4, g u +}^{k, N O D} (f_{0, g u -}^{k + 1} {\bar{ϕ}}_{0, g}^{k + 1} + f_{0, g u +}^{k} {\bar{ϕ}}_{0, g}^{k}) \end{array}

(13b)

where $D_{2, g u +}^{k, FDM}$ and $D_{2, g u +}^{k, FDM}$ are the zeroth-order and second-order pseudo-diffusion coefficients, respectively; $D_{2, g u +}^{k, FDM}$ and $D_{2, g u +}^{k, FDM}$ are the zeroth-order and second-order coupling corrective coefficients, respectively; and $f_{0, g u -}^{k + 1}$ and $f_{0, g u -}^{k + 1}$ are the zeroth-order and second-order discontinuity factors, respectively.

A new nonlinear iterative method [25] was proposed based on the above mentioned coupling corrective relationship, and the NLSP3 reactor core code was developed. Figure 2 illustrates the nonlinear iterative strategy. Using the NLSP3 code and global cross sections, the global flux distribution can be obtained easily.

Fig. 2

Nonlinear iterative strategy of the NLSP3 code

2.3

GVR for MC Criticality Calculation

The MC criticality calculation is composed of several inactive and active cycles. Each cycle is a fixed-source calculation with a cycle-estimated result. All the active cycles are averaged as follows to obtain the final result of the MC criticality calculation:

\bar{X} = \frac{\sum_{n = 1}^{N} \bar{Y_{n}}}{N},

(14)

where the subscript n denotes the cycle index, $\bar{X}$ the final result of the MC criticality calculation, N the total number of active cycles, and $\bar{Y_{n}}$ the cycle-estimated result of the n^th active cycle. $\bar{Y_{n}}$ is given by

\bar{Y_{n}} = \frac{\sum_{i = 1}^{M} Y_{n}^{i}}{M},

(15)

where $Y_{n}^{i}$ is the i^th neutron history of the n^th active cycle, and M is the number of fission sources in a single cycle.

The standard deviation of the final result is calculated from the deviation between the individual active cycles:

S_{\bar{X}} = \sqrt{\frac{S^{2}}{N}} = \sqrt{\frac{\bar{X^{2}} - {\bar{X}}^{2}}{N - 1}},

(16)

where

\bar{X^{2}} = \frac{{\sum_{n = 1}^{N} \bar{Y_{n}}}^{2}}{N} .

(17)

Therefore, the standard deviation increases with the deviation between the individual active cycles. The deviation between the individual active cycles should therefore be reduced to reduce the standard deviation of the MC criticality calculation. As stated before, a fixed-source calculation is performed for each cycle, and the deviation between the individual cycles is reduced by various VR techniques.

Suppose that all the active cycles are performed with the same weight windows. The individual cycle-estimated results will then be close to another. This leads to a reduction in the standard deviation. Therefore, the key to improving the MC criticality calculation is to obtain high-precision weight window parameters efficiently.

In the SP3-coupled GVR method, only the forward SP3 calculation is performed. Because there is no known adjoint theory of the SP3 method, we cannot obtain the adjoint flux for the time being. We have therefore focused on GVR for the criticality calculations in the meantime. We believe that the forward flux distribution is sufficient for determining the weight window parameters in place of the fixed-source calculation. Therefore, we use the forward flux distribution to approximate the weight window parameters.

The global flux distribution determined by the SP3 equation is denoted as $ϕ_{g}^{k}$ , where the superscript k denotes the spatial region and the subscript g denotes the energy interval. The SP3-coupled method can only bias the spatial position of the neutrons. Therefore, the independent energy variable should be eliminated by accumulating the energy-dependent flux in the spatial region:

ϕ_{k} = \sum_{g = 1}^{G} ϕ_{g}^{k},

(18)

where φ_k is the space-dependent flux distribution, which is normalized before it is used to approximate the global weight windows, and

{\bar{ϕ}}_{k} = ϕ_{k} / max (ϕ_{k}),

(19)

where max(φ_k) is the largest space-dependent flux. Using the normalized global flux distribution, the global weight windows are determined as

W_{s, k} = {\bar{ϕ}}_{k},

(20)

where W_s,k is the global space-dependent survival weight window in the k^th space region. Using Eq. (20), the weight windows of all the spatial regions can be determined, following which the ceiling and floor weight windows can be determined using the following formulas:

W_{c, k} = \frac{5 W_{s, k}}{3},

(21)

and

W_{f, k} = \frac{W_{s, k}}{3} .

(22)

where W_c,k and W_f,k are the ceiling and floor weight windows in the k^th space region, respectively. The standard weight window with the weight windows of all the spatial regions is used for GVR.

NUMERICAL RESULTS

The SP3-coupled GVR method was tested on the C5G7 benchmark [28]. In the original C5G7 model, GVR can be easily achieved because the water shield layer is not sufficiently thick to pose a deep penetration problem. To better demonstrate the advantages of the SP3-coupled method for GVR, we artificially increased the thickness of the shield layer. The materials and geometry of the reactor core were derived from the C5G7 benchmark, and the core region was surrounded by a water shield layer, which was 64.362 cm thick and three times thicker than the side of the core assembly side as shown in Fig. 3.

Fig. 3

(Color online) Modified C5G7 benchmark

The superimposed mesh geometry was used for SP3-coupled GVR simulation. The C5G7 reactor core had a 4×4 assembly layout and the assemblies was surrounded by a water shield that was three times thicker as the assembly. The model was divided into a total of 10×10=100 material mesh in the modified C5G7 benchmark. The MGCS of the 100 mesh were determined through the MC global homogenization over the four energy intervals of ［0, 6.25×10^-7］, [6.25×10^-7, 5.53×10^-3], [5.53×10^-3, 8.21×10^-1], [8.21×10^-1, 20] MeV [29].

A fixed-source calculation was performed for the MC global homogenization. The choice of the number of source neutrons per material mesh is a challenging issue. A greater number of source neutrons for the fixed-source calculation will increase the precision of the obtained cross sections but will also require more computation. To achieve a balance between efficiency and accuracy, the number of source neutrons per material mesh was chosen as 10,000 based on our experience.

Therefore, a total of 1,000,000 source neutrons were simulated for the MC global homogenization. The SP3 solution only provided the weight windows, and the final result was given by the MC criticality calculation. Thus, even if there was an error in the MGCS results, the accuracy of the final result was not affected. This is why the number of source neutrons per material mesh was empirically determined for the solution of the MGCS.

The RMC code called the SP3 solver for the global calculation with the cross-sections of the 10×10 material mesh. Each material mesh was divided into a total of 4×4=16 nodes with each node being 5 cm × 5 cm in size. Therefore, a total of 40×40=1600 nodes were used to solve the SP3 equation. After solving the SP3 equation, the flux distribution in the 1600 nodes was obtained. The three convergence criteria for the inner, outer, and k_eff iterations in the NLSP3 code were set to 10^-3, 10^-3, and 10^-4, respectively. The calculations were performed with 24 threads, and the calculation time taken for the MC global homogenization and SP3 solver was 1.14 minutes.

Given the global flux distribution, the global weight windows can be determined using Eqs. (20) – (22). The MC criticality calculation was then performed using the global weight windows. There were 10,000 neutrons per cycle, 50 inactive cycles, and 450 active cycles in the SP3-coupled GVR criticality calculation. 15.67 minutes were required to complete the MC global homogenization, SP3 solver solution, and GVR criticality calculation.

For comparison, the standard MC simulation method was also performed to calculate the flux distribution of the 1600 nodes. Because GVR was not performed in the standard simulation method, we call this the direct simulation. There were 1000,000 neutrons per cycle, 50 inactive cycles, and 450 active cycles in the simulation. The calculation time was 151.80 minutes. The distributions of the relative errors in the SP3-coupled GVR and the direct simulations are compared in Fig. 4. The distributions of the total flux in the SP3-coupled GVR and the direct simulations are compared in Figure 5.

Fig. 4

(Color online) Variance distribution of the SP3-coupled GVR and direct simulations

Fig. 5

(Color online) Distribution of the total flux in the SP3-coupled GVR and direct simulations

As shown in Fig. 5, the distributions of the total flux obtained by the direct and the SP3-coupled GVR simulations are in good agreement. The deviation was calculated to be less than three times the standard deviation. Meanwhile, it can be seen from Figure 4 that when the modified C5G7 benchmark was simulated directly, the relative errors increased with the distance away from the reactor core. Compared with the direct simulation, the global variance distribution in the SP3-coupled GVR simulation was flattened in the entire system. High accuracy was achieved in both the reactor core and the shielding layer. Therefore, the SP3-coupled GVR method can improve the efficiency of MC criticality calculations. To better quantify its GVR performance, three quantities were compared:

1. The average figure of merit (AV.FOM), which is defined as

A V . F O M = \frac{N}{T \cdot \sum_{i = 1}^{N} {Re}_{i}^{2}} .

(23)

2. The standard deviation of the relative deviation, which is defined as

σ_{Re} = \sqrt{\frac{1}{N} \cdot \sum_{i = 1}^{N} {Re}_{i}^{2} - \frac{1}{N^{2}} \cdot {(\sum_{i = 1}^{N} {Re}_{i}^{2})}^{2}} .

(24)

3. The ratio of the maximum to the minimum variances (σ_max/σ_min), i.e.,

σ_{\max} / σ_{\min} = \frac{\max ({Re}_{i})}{\min ({Re}_{i})} .

(25)

The quantitative results for the direct and SP3-coupled GVR simulations are compared in Table 1.

Results of GVR and direct simulations

Quantities	Direct simulation	SP3-coupled GVR
AV.FOM	1.1364×10¹	3.8394×10²
σ.Re	2.4070×10^-2	1.2893×10^-2
σ_max/σ_min	6.1998×10²	1.2196×10¹
Times (min)	151.80	15.67

The results show that, compared with the direct simulation, the average figure of merit (AV. FOM) was improved 33 folds in the SP3-coupled GVR method. Meanwhile, the ratio of the maximum to minimum variances was reduced by a factor of five. The value of the SP3-coupled GVR method for GVR and its ability to improve the efficiency of the MC criticality calculation are hence clearly demonstrated.

It has been argued that it is unnecessary to couple the SP3 method with the MC method for GVR, and that the diffusion method is sufficient. To demonstrate the advantage of coupling the SP3 method with the MC method for GVR over coupling the diffusion method, we developed a diffusion-coupled GVR method based on the RMC code. The flowcharts of the diffusion-coupled GVR and the SP3-coupled GVR methods are almost identical except for the solver used to solve the global flux distribution. In the SP3-coupled GVR method, the SP3 equation is solved, while in the diffusion-coupled method, the diffusion equation is solved. To solve the diffusion equation, the diffusion coefficient is required. Therefore, the total cross-section $Σ_{t, g}$ used in the SP3 equation was converted into the diffusion coefficient D_g.

Calculations were then performed on the same model using the diffusion-coupled GVR method with the same calculation parameters. The obtained AV.FOM from the results of the diffusion-coupled method was 136.81. As shown in Table 1, the AV.FOM of the result of the SP3-coupled method was 383.94, which is 2.81 times that of the diffusion-coupled method. Therefore, the SP3-coupled GVR method showed better GVR performance than the diffusion-coupled GVR method. We believe that the SP3 method has great potential to be coupled with the MC method for accelerated calculations and GVR.

DISCUSSION AND CONCLUSION

An automatic process based on MC global homogenization, the SP3 global solver, and GVR criticality calculation was developed in the SP3-coupled GVR method. The numerical results show that this method can improve the efficiency of the MC criticality calculation.

The SP3-coupled method proposed in this paper was developed for MC criticality calculations. It therefore superficially resembles the traditional two-step reactor physical calculation framework. However, the SP3-coupled method differs from the traditional framework in three aspects: (1) global homogenization is performed, instead of assembly homogenization; (2) the homogenization process is completed during the fixed-source calculation, and not during the criticality calculation; (3) global MC homogenization and the global reactor calculation are achieved in a one-step process instead of a two-step process. Moreover, only the weight window parameters for the MC criticality calculation are output in the SP3-coupled method. The final results are obtained from the MC calculation rather than the global calculation of the NLSP3 code.

The SP3-coupled method has similarities and differences with other hybrid MC methods, such as the fission matrix and the response matrix methods. They share similar calculation processes in which some local information of the system is recorded in advance and used to build a global formula, and the global information is obtained by solving the global formula. In the fission matrix method, the fission coupling relationship between different regions is recorded first. In the response matrix method, the relationship between the currents and sources in different regions is recorded first. In the SP3-coupled method, the homogenized cross-sections in different regions are recorded first. The matrix relationship is solved in the fission matrix and response matrix methods to obtain global information. In contrast, in the SP3-coupled method, the neutron transport equation is solved. We therefore believe that the SP3-coupled method has advantages in terms of numerical stability, development potential, and application scope.

In this study, the traditional SP3 equation is solved in the SP3-coupled method and the forward flux is used to approximate the weight windows. To further improve its performance for GVR, the rigorous SPN theory proposed by Prof. Yung-An Chao and the adjoint method for the SPN theory should be studied to obtain better weight window parameters.

Because the SP3-coupled method is helpful for the GVR of the MC criticality calculation, we hope to apply it to reactor shielding calculations in the future. We are attempting to perform reactor shielding design using MC criticality and burnup calculations, and hope to have an accurate radiation source for reactor shielding calculations to analyze the influence of the reactor dynamics on the reactor shielding.

REFERENCE

[1]

Ma Y., Liu S., Luo Z. et al.,

RMC/CTF multiphysics solutions to VERA core physics benchmark problem 9

. Ann. Nucl. Energy 133, 837-852 (2019). doi: 10.1016/j.anucene.2019.07.033.