The batch method

The batch method de-correlates the tally results by combining a certain number of active cycles into batches. In the MC simulation, $Q$ is assumed to be the physical quantity of interest. Next, its conventional MC estimate value $\overline{Q}$ and estimate variance ${\sigma }_s^2\left(\overline{Q}\right)$ are assumed to be [3] $\begin{array}{c}\overline{Q}=\frac{1}{N}{\displaystyle \sum _{i=1}^N{Q}^{\left(i\right)}}\end{array}$ (1) $\begin{array}{c}{\sigma }_{\text{s}}^2\left(\overline{Q}\right)=\frac{1}{N\left(N-1\right)}{\displaystyle \sum _{i=1}^N{\left(Q^{\left(i\right)}-\overline{Q}\right)}^2}\end{array}$ (2) where $Q^{\left(i\right)}$ is the estimate in cycle i, and $N$ is the number of active cycles.

The apparent variance of $\overline{Q}$, ${\sigma }_{\text{a}}^2\left(\overline{Q}\right)$ is defined as the mathematical expectation of the estimated variance. $\begin{array}{c}{\sigma }_{\text{a}}^2\left(\overline{Q}\right)=E\left({\sigma }_{\text{s}}^2\left(\overline{Q}\right)\right)=\frac{1}{N}E\left({\left(Q^{\left(i\right)}\right)}^2\right)-\frac{1}{N}{E}^2\left(Q^{\left(i\right)}\right)-\frac{1}{N^2}\frac{1}{N-1}{\displaystyle \sum _i{\displaystyle \sum _{j\ne i}\text{cov}}}\left(Q^{\left(i\right)},Q^{\left(j\right)}\right)\end{array}$ (3)

If we define ${\sigma }^2\left(Q^{\left(i\right)}\right)=E\left({\left(Q^{\left(i\right)}\right)}^2\right)-E^2\left(Q^{\left(i\right)}\right)$, then Eq. (3) can be simplified as $\begin{array}{c}{\sigma }_{\text{a}}^2\left(\overline{Q}\right)=\frac{1}{N}{\sigma }^2\left(Q^{\left(i\right)}\right)-\frac{1}{N^2}\frac{1}{N-1}{\displaystyle \sum _i{\displaystyle \sum _{j\ne i}\text{cov}}}\left(Q^{\left(i\right)},Q^{\left(j\right)}\right)\end{array}$ (4)

The true variance ${\sigma }_{\text{r}}^2\left(\overline{Q}\right)$ of $\overline{Q}$ is expressed by Eq. (5). $\begin{array}{c}{\sigma }_{\text{r}}^2\left(\overline{Q}\right)=E\left({\overline{Q}}^2\right)-{\left(E\left(\overline{Q}\right)\right)}^2=\frac{1}{N}{\sigma }^2\left(Q^{\left(i\right)}\right)+\frac{1}{N^2}{\displaystyle \sum _i{\displaystyle \sum _{j\ne i}\text{cov}}}\left(Q^{\left(i\right)},Q^{\left(j\right)}\right)\end{array}$ (5)

The difference between ${\sigma }_{\text{r}}^2\left(\overline{Q}\right)$ and ${\sigma }_{\text{a}}^2\left(\overline{Q}\right)$ is $\begin{array}{c}{\sigma }_{\text{r}}^2\left(\overline{Q}\right)-{\sigma }_{\text{a}}^2\left(\overline{Q}\right)=\frac{1}{N\left(N-1\right)}{\displaystyle \sum _i{\displaystyle \sum _{j\ne i}\text{cov}}}\left(Q^{\left(i\right)},Q^{\left(j\right)}\right).\end{array}$ (6)

Thus, the correlation between the tally results of every cycle is the direct cause of the variance underestimation problem.

Furthermore, we used the batch method by combining a certain number of active cycles into batches and tallying the physical quantities at the end of each batch. Suppose that the batch length is $l_{\text{b}}$ and the number of batches is $n_{\text{b}}$. The coefficients confront Eq. (7). $\begin{array}{c}{l}_{\text{b}}\times n_{\text{b}}=N \end{array}$ (7)

We defined $Q_{\text{b}}^{\left(j\right)}$ as the estimate in batch j and $Q_{\text{b}}^{j\left(i\right)}$ as the estimate of cycle i in batch j. The estimated value ${\overline{Q}}_{\text{b}}$ and sample variance ${\sigma }_{\text{bs}}^2\left({\overline{Q}}_{\text{b}}\right)$ can be as Eq. (8-9). $\begin{array}{c}{\overline{Q}}_{\text{b}}=\frac{1}{n_{\text{b}}}{\displaystyle \sum _{j=1}^{n_{\text{b}}}{Q}_{\text{b}}^{\left(j\right)}}=\frac{1}{n_{\text{b}}}{\displaystyle \sum _{j=1}^{n_{\text{b}}}\left(\frac{1}{l_{\text{b}}}{\displaystyle \sum _{i=1}^{l_{\text{b}}}{Q}_{\text{b}}^{j\left(i\right)}}\right)}=\frac{1}{N}{\displaystyle \sum _{i=1}^N{Q}^{\left(i\right)}}=\overline{Q}\end{array}$ (8) $\begin{array}{c}{\sigma }_{\text{bs}}^2\left({\overline{Q}}_b\right)=\frac{1}{n_{\text{b}}\left(n_{\text{b}}-1\right)}{\displaystyle \sum _{i=1}^{n_{\text{b}}}{\left(Q_{\text{b}}^{\left(i\right)}-\overline{Q}\right)}^2}\end{array}$ (9)

The apparent variance of ${\overline{Q}}_{\text{b}}$, ${\sigma }_{\text{ba}}^2\left(\overline{Q}\right)$ is defined as the expected value of the sample variance. $\begin{array}{c}{\sigma }_{\text{ba}}^2\left({\overline{Q}}_b\right)=E\left({\sigma }_{\text{bs}}^2\left({\overline{Q}}_{\text{b}}\right)\right)=\frac{1}{n_{\text{b}}}{\sigma }^2\left(Q_b^{\left(i\right)}\right)+\frac{1}{n_{\text{b}}-1}\frac{1}{n_{\text{b}}^2}{\displaystyle \sum _i{\displaystyle \sum _{j\ne i}\text{cov}}}\left(Q_{\text{b}}^{\left(i\right)},Q_{\text{b}}^{\left(j\right)}\right)\end{array}$ (10) where ${\sigma }^2\left(Q_{\text{b}}^{\left(i\right)}\right)=E\left({\left(Q_{\text{b}}^{\left(i\right)}\right)}^2\right)-E^2\left(Q_{\text{b}}^{\left(i\right)}\right)$.

Then, the difference between ${\sigma }_{\text{br}}^2\left(\overline{Q}\right)$ and ${\sigma }_{\text{ba}}^2\left(\overline{Q}\right)$ is $\begin{array}{c}{\sigma }_{\text{br}}^2\left(\overline{Q}\right)-{\sigma }_{\text{ba}}^2\left(\overline{Q}\right)=\frac{1}{n_{\text{b}}\left(n_{\text{b}}-1\right)}{\displaystyle \sum _i{\displaystyle \sum _{j\ne i}\text{cov}}}\left(Q_{\text{b}}^{\left(i\right)},Q_{\text{b}}^{\left(j\right)}\right)\end{array}$ (11)

The correlation among batch tally results is generally less significant than that among cycle tally results in the PI process, owing to the weaker FSD inter-batch correlation. Thus, the difference between ${\sigma }_{\text{br}}^2\left(\overline{Q}\right)$ and ${\sigma }_{\text{ba}}^2\left(\overline{Q}\right)$ was smaller than that between ${\sigma }_{\text{r}}^2\left(\overline{Q}\right)$ and ${\sigma }_{\text{a}}^2\left(\overline{Q}\right)$. However, setting an appropriate batch length is challenging. A large batch length can almost eliminate underestimation, whereas the computational cost can be too high because it requires a minimum number of batches in common practice.

A small batch length is appropriate for practical use, but the magnitude of variance underestimation elimination is unknown. Determining the optimum batch length was of interest in this study.

Sliced Wasserstein distance algorithm

The Sliced Wasserstein (SW) distance was developed from Optimum Transport theory [21]. Its initial definition is the minimum "cost" to transform the 1D distribution $p$ into $q$, as shown in Fig. 1. It can be efficiently calculated by ${\text{SW}}_k\left(p,q\right)=\frac{1}{N}{\displaystyle \sum _{i=1}^N\text{d}}\left(x_i,y_i\right),$ (12) where k is the X/Y/Z direction, $p,q$ are two 1D distributions, $x_i\in p$, $y_i\in q$ are the coordinates after sorting the neutron points in the X/Y/Z direction, and N is the number of particles per cycle in the MC simulations.

Fig. 1Full-screen View

Minimum "cost" to transform distribution $p$ into $q$

In this paper, the 1-norm distance is used to calculate $d\left(x,y\right)$: $\begin{array}{c}\text{d}\left(x,y\right)=\Vert x-y_1\Vert =\sum \left|x_i-y_i\right|\end{array}.$ (13)

The SW distance algorithm is a commonly used mesh-free algorithm for measuring the difference between discrete vectors. Guo et al. (2021) introduced this algorithm into the scope of MC source convergence diagnosis [22]. Shen et al. (2022) adopted this algorithm to measure fission source error in the MC PI process [23]. In the calculation, the sorting procedure uses the fast sorting algorithm costs $O\left(N\text{log}\left(N\right)\right)$ operations, and the SW distance calculation costs $O\left(N\right)$ operations. Thus, the total cost is $O\left(N\text{log}\left(N\right)\right)$. The efficiency penalty is approximately 2% in MC criticality calculations for the PWR full-core model [22].

SW-based FSD inter-cycle correlation observation method

The PI method in the MC algorithm is expressed by Eq. (14) in the discretization. $\begin{array}{c}{s}^{\left(i+1\right)}=\frac{H{s}^{\left(i\right)}}{k^{\left(i\right)}}+{ϵ}^{\left(i\right)}\end{array},$ (14) where $s^{\left(i\right)}$ is the fission source vector in cycle I, $H$ is the fission matrix, ${ϵ}^{\left(i\right)}$ is the stochastic error component, and $k^{\left(i\right)}$ is the estimation of fundamental mode eigenvalue in cycle i.

The number of particles per cycle is defined as m, the fundamental mode eigenvalue, and the normalized eigenvector of fission matrix $H$ is defined as $k_0$ and $s_0$. Eq. (15) can be expressed as $\begin{array}{c}H{s}_0=k_0{s}_0\end{array}.$ (15)

The fission source error in cycle i is introduced as $\begin{array}{c}{e}^{\left(i\right)}=s^{\left(i\right)}-m{s}_0.\end{array}$ (16)

By substituting and considering the magnitude of the error terms [10], we disregarded some terms and obtained $\begin{array}{c}{e}^{\left(i\right)}\cong A{e}^{\left(i-1\right)}-\frac{\tau H}{k_{0}m}{e}^{\left(i-1\right)}A{e}^{\left(i-1\right)}+{ϵ}^{\left(i-1\right)} ,\end{array}$ (17) where $\tau$ is (1, 1, …, 1), and $A$ is the noise propagation matrix, written as Eq. (18). $\begin{array}{c}A=\frac{I-s_0\tau }{k_0}H\end{array}$ (18)

We iterated Eq. (17) i times and obtained Eq. (19). $\begin{array}{c}{e}^{\left(i\right)}\cong A^i{e}^{\left(0\right)}+O\left(m^{-1}\right)+{\displaystyle \sum }_{j=1}^{i-1}{A}^j{ϵ}^{\left(i-j\right)} \end{array}$ (19)

The fission source error consists of three components: the initial error, bias error, and random error terms. The FSD random error term is the dominant component after the source convergence, as expressed in Eq. (20) [24]. $\begin{array}{c}{e}_{\text{R}}^{\left(i\right)}\cong {\displaystyle \sum _{j=0}^{i-1}{A}^j{ϵ}^{\left(i-j\right)}} \end{array},$ (20) where $e_R^{\left(i\right)}$ is the fission source random error term of cycle i, $A$ is the error propagation matrix, and ${ϵ}^{\left(j\right)}$ is the stochastic error component induced in cycle j.

The FSD random error term is inter-cycle correlated because it is composed of stochastic error components induced in adjacent cycles. This is the theoretical basis of the FSD inter-cycle correlation.

The fission error difference between adjacent cycles can be written as $\begin{array}{c}{e}_{\text{R}}^{\left(i\right)}-e_{\text{R}}^{\left(i-1\right)}={ϵ}^{\left(i\right)}+{\displaystyle \sum _{j=1}^{i-1}\left(A-I\right)}{A}^{j-1}{ϵ}^{\left(i-j\right)},\end{array}$ (21) where $I$ is the identity matrix.

In fission systems, noise propagation matrix $A$ is similar to the identity matrix. By comparing the scalar error differences, we obtained $\begin{array}{c}\Vert e_{\text{R}}^{\left(i\right)}\Vert >\Vert e_{\text{R}}^{\left(i\right)}-e_{\text{R}}^{\left(i-1\right)}\Vert \end{array},$ (22) where ${\epsilon }^{\left(i\right)}$ is the scalar fission source error in cycle i.

Furthermore, the fission error difference between cycles can be written as $\begin{array}{c}{e}_{\text{R}}^{\left(i\right)}-e_{\text{R}}^{\left(i-j\right)}={\displaystyle \sum _{k=0}^{j-1}{A}^k{ϵ}^{\left(i-k\right)}}+{\displaystyle \sum _{t=j}^{i-1}\left(A^t-A^{t-j}\right)}{ϵ}^{\left(i-t\right)}.\end{array}$ (23)

Similarly, the scalar error differences are expected to be $\Vert e_{\text{R}}^{\left(i\right)}\Vert \begin{array}{c}>\Vert e_{\text{R}}^{\left(i\right)}-e_{\text{R}}^{\left(i-j\right)}\Vert >\dots >\Vert e_{\text{R}}^{\left(i\right)}-e_{\text{R}}^{\left(i-1\right)}\Vert ,j>1\end{array}$ (24)

The scalar error differences increase with the interval j. When $E\left(\Vert e_R^{\left(i\right)}\Vert \right)\approx E\left(\Vert e_R^{\left(i\right)}-e_R^{\left(i-n\right)}\Vert \right)$, the FSD inter-cycle correlation can be ignored after n cycles.

The SW-based FSD inter-cycle correlation observation method uses the combined Sliced Wasserstein distance in three directions to measure the scalar fission source error ${\epsilon }^{\left(i\right)}$, as expressed in Eq. (25). $\begin{array}{c}\text{SW}\left(s^{\left(i\right)},s^{\left(j\right)}\right)\text{ is used to measure }\Vert e_R^{\left(i\right)}-e_R^{\left(j\right)}\Vert ,\end{array}$ (25) where $s^{\left(i\right)}$ is the fission source distribution in cycle i.

The combined SW distance is calculated in Eq. (26). $\begin{array}{c}\text{SW}\left(s^{\left(i\right)},s^{\left(j\right)}\right)={\text{SW}}_x\left(s^{\left(i\right)},s^{\left(j\right)}\right)\times {\text{SW}}_y\left(s^{\left(i\right)},s^{\left(j\right)}\right)\times {\text{SW}}_z\left(s^{\left(i\right)},s^{\left(j\right)}\right)\end{array}$ (26)

We combined Eqs. 24 and 25, and the average SW distance between the FSD of adjacent cycles increased along with the adjacent cycle number, as expressed in Eq. (27). $\begin{array}{c}\forall i>k>j,E\left(\text{SW}\left(s^{\left(i\right)},s^{\left(j\right)}\right)\right)\ge E\left(\text{SW}\left(s^{\left(k\right)},s^{\left(j\right)}\right)\right)\end{array}$ (27)

The FSD inter-cycle correlation can be ignored when the average FSD SW distance fluctuates.

Inter-cycle correlation length and Optimum batch length

We defined the interval cycle length $n_0$ where the average FSD SW distance starts to fluctuate as the FSD inter-cycle correlation length. The interval length $n_0$ conforms to Eq. (28). $\begin{array}{c}\forall \left( j\ge n_0 \right), E\left(\overline{\text{SW}\left(s^{\left(i\right)},s^{\left(i-n_0\right)}\right)}\right)\cong E\left(\overline{\text{SW}\left(s^{\left(i\right)},s^{\left(i-j\right)}\right)}\right)\end{array}$ (28)

According to Eq. (25), the FSD SW distance was used to measure the scalar fission source error. Thus, Eq. (29) can be expressed as $\begin{array}{c}\forall \left( j\ge n_0 \right),\text{} E\left(\overline{e^{\left(i\right)}-e^{\left(i-n_0\right)}}\right)\cong E\left(\overline{e^{\left(i\right)}-e^{\left(i-j\right)}}\right).\end{array}$ (29)

If the interval cycle length j is greater than $n_0$, the correlation between the FSDs of cycle i and cycle $i-j$ declines statistically to negligible.

According to the Law of Large Numbers, we can obtain $\begin{array}{c}\forall \left(j\ge n_0\right),E\left({\displaystyle \sum _{k=0}^l\Vert e^{\left(i+k\right)}-e^{\left(i+k-n_0\right)}\Vert }\right)\cong E\left({\displaystyle \sum _{k=0}^l\Vert e^{\left(i+k\right)}-e^{\left(i+k-j\right)}\Vert }\right),\end{array}$ (30) where l is group length.

We assumed the cumulated scalar fission error of l cycles to be uncorrelated after $n_0$ cycles. Specifically, when group length l is set as the FSD inter-cycle correlation length $n_0$, Eq. (31) can be expressed as $\begin{array}{c}\forall \left(j\ge n_0\right),E\left({\displaystyle \sum _{k=0}^{n_0}\Vert e^{\left(i+k\right)}-e^{\left(i+k-n_0\right)}\Vert }\right)\cong E\left({\displaystyle \sum _{k=0}^{n_0}\Vert e^{\left(i+k\right)}-e^{\left(i+k-j\right)}\Vert }\right).\end{array}$ (31)

The cumulative scalar fission error of $n_0$ cycles is statistically uncorrelated after $n_0$ cycles. The flux results $\left\{Q_{\text{b}}^{\left(i\right)}\right\}$ are then statistically uncorrelated using $n_0$ as the batch length. The variance underestimation problem can be solved consequently, as expressed in Eq. (32). $\begin{array}{c}{\sigma }_{\text{br}}^2\left(\overline{Q}\right)-{\sigma }_{\text{ba}}^2\left(\overline{Q}\right)=\frac{1}{n_{\text{b}}\left(n_{\text{b}}-1\right)}{\displaystyle \sum _i{\displaystyle \sum _{j\ne i}\text{cov}}}\left(Q_{\text{b}}^{\left(i\right)},Q_{\text{b}}^{\left(j\right)}\right)\cong 0\end{array}$ (32)

In conclusion, the FSD inter-cycle correlation length is the optimum batch length.

Workflow

The SeVUE method comprises two procedures: the Estimation procedure and the Elimination procedure.

In the Estimation procedure, the SeVUE method requires an MC calculation to estimate the FSD inter-cycle correlation length of each model. This calculation should be in the criticality mode with no tallies. The number of particles per cycle can be set to a relatively small value because the FSD inter-cycle correlation is almost independent of it. For example, setting 10,000 neutrons per cycle is sufficient for a full-core model. By contrast, the number of active cycles should be relatively large. Thousands of active cycles are recommended for models with a high dominance ratio. The SW distances between FSDs can be calculated using external scripts (Fig. 4) or embedded codes. At the end of this procedure, the relationship between the average FSD SW distances and the interval cycle length can be plotted (Fig. 3). The FSD inter-cycle correlation length is the point at which the average FSD SW distance starts to fluctuate.

Fig. 4Full-screen View

Coupled workflow to estimate FSD inter-cycle correlation length

Fig. 3Full-screen View

Relationship between FSD average SW distances and interval cycle length

In the Elimination procedure, the batch method was used to eliminate variance underestimation. According to Eq. (19-22), the FSD inter-cycle correlation length was set as the batch length. The number of active cycles should be at least 10 times the FSD inter-cycle correlation length. The calculated variances of the tally results were close to the actual variances.

Implementation

The SeVUE method was implemented using the MC code RMC [20]. RMC is a neutron photon-electron transport code developed by the Reactor Engineering Analysis Lab (REAL) of Tsinghua University. RMC has been applied for nuclear reactor criticality, burnup, shielding, and neutronics-thermal coupling calculations [25-28].

In the Estimation procedure, an external script was used to calculate the SW distances between FSDs of different cycles. In the RMC code, the fission source distribution is the output for each cycle. The external script can run independently and be easily coupled with other MC codes. The coupled workflow is illustrated in Fig. 3.

In the Elimination procedure, the batch method is implemented in the RMC code. Tally results were processed at the end of each batch.

Flux distribution of sphere array model

The sphere array model was obtained from the OECD/NEA source convergence benchmark, with a dominance ratio of approximately 0.886. It comprised of 25 fissile metal spheres in a $5\times 5\times 1$ array (Fig. 5). The metal spheres were filled with highly enriched uranium. The center sphere was larger than the other eight spheres, with a radius of 10 cm. The flux distribution in this model was of interest in this study. The number of particles per cycle of 5,000, 50,000, and 5,000,000 was analyzed to observe the FSD inter-cycle correlation lengths. Each calculation had 1,000 inactive and 20,000 active cycles.

Fig. 5Full-screen View

(Color online) Geometric configuration of sphere array model

4.1.1

Reference flux distribution and variance underestimation phenomenon

The reference flux distribution and relative standard deviation (RSD) distributions were obtained using 100 independent calculations. The fluxes in each sphere were tallied. The reference flux results and RSDs are shown in Fig. 6. According to the Law of large numbers, the RSD of the sample variance was approximately $\sqrt{\frac{2}{N-1}}$ (where N is the number of samples). Thus, the reference standard deviation has a $\pm 10\text{\%}$ fluctuation.

Fig. 6Full-screen View

(Color online) The reference flux results and RSDs of the Sphere array model. (a) Flux results; (b) Relative standard deviation.

The variance underestimation ratio distribution in Fig. 7 was obtained by comparing the relative standard deviations in a single simulation with the reference results. The Variance Underestimation Ratio (VUR) is defined in Eq. (33). $\begin{array}{c}\text{VUR}=\frac{{\sigma }_{\text{r}}-{\sigma }_{\text{a}}}{{\sigma }_{\text{r}}}\times 100\%,\end{array}$ (33) where ${\sigma }_a$ is the calculated standard deviation, and ${\sigma }_r$ is the true standard deviation.

Fig. 7Full-screen View

(Color online) RSDs and VURs distribution of a single MC simulation in sphere array model. (a) Relative standard deviation; (b) VURs

As shown in Fig. 7, the VURs in each sphere varied from 24% to 47%. The average VUR was 37.3%.

4.1.2

Estimation procedure

The number of particles per cycle of 5,000, 50,000, and 5,000,000 was analyzed to observe the FSD inter-cycle correlation lengths in this procedure. Fig. 8 shows that the FSD average SW distance evolves with the interval cycle length. A visual check of the curves showed that the FSD inter-cycle correlation lengths were approximately 50 for all three batch sizes. Thus, the FSD inter-cycle correlation was almost independent of batch size.

Fig. 8Full-screen View

FSD average SW distance evolves with interval cycle length in sphere array model. (a) 5,000 neutrons per cycle; (b) 50,000 neutrons per cycle; (c) 500,000 neutrons per cycle.

Semi-quantitative information can be obtained from Fig. 8. The initial FSD average SW distance was approximately 40% of the stationary SW distance. This "similarity coefficient" can be further used to measure the variance underestimation ratio of the total flux in this model.

4.1.3

Elimination procedure

We verified Eqs. (28-32) by setting different batch lengths. The batch lengths are listed in Table. 1. Notably, in the SeVUE method, only the case with $l_{\text{b}}=50$ should be calculated.

Table 1Full-screen View

Parameters of batch method in sphere array model

Serial number	Batch length/cycle	Number of batches	Number of active cycles
1	1	20,000	20,000
2	2	10,000	20,000
3	4	5000	20,000
4	5	4000	20,000
5	8	2500	20,000
6	10	2000	20,000
7	16	1250	20,000
8	20	1000	20,000
9	25	800	20,000
10	32	625	20,000
11	40	500	20,000
12	50	400	20,000

Show more

Figure 9 shows that the average VURs evolve with the batch length. The average VURs of different batch sizes are approximately 37% when the batch length is one cycle. With the batch length increases, the average VURs gradually declines. When batch length is 50, the average VURs are approximately $\pm 1\%$. Thus, 50 is the proper batch length for this model. This value is in satisfactory agreement with the FSD inter-cycle correlation length estimated in the former procedure.

Fig. 9Full-screen View

Average VUR evolves with batch length in sphere array model

The results of a single calculation for $l_{\text{b}}=50$ are shown in Fig. 10. The VUR of each sphere varied from -5.0% to 2.6%, with an average VUR of -1.9%. Considering the fluctuations in the reference standard deviations, these differences were acceptable.

Fig. 10Full-screen View

(Color online) Relative standard deviation and VURs distributions of an MC simulation with $l_{\text{b}}=50$. (a) Relative standard deviation; (b) VURs.

Power distribution of BEAVRS 3D full-core model

The BEAVRS full-core model is a detailed reactor core model used to test indicators in real-world PWRs. The geometric configuration is illustrated in Fig. 11. It consists of 193 fuel assemblies in a 15 $\times$ 15 $\times$ 1 array with a dominance ratio of 0.989. The power distribution in this model was of interest in this study. The number of particles per 50,000 cycles was analyzed to observe the FSD inter-cycle correlation length. Each calculation had 1,000 inactive and 20,000 active cycles. The total number of neutron histories in the active cycles was 1 billion, which is practical for actual calculations.

Fig. 11Full-screen View

(Color online) Geometric configuration of 3D BEAVRS full-core model. (a) Front view; (b) Left view.

4.2.1

Reference power distribution and variance underestimation phenomenon

The reference power distribution and standard deviation distribution were obtained using 30 independent calculations, similar to the sphere array model. The power of each assembly was tallied. The reference results and RSDs are shown in Fig. 12.

Fig. 12Full-screen View

(Color online) Reference power results and relative standard deviations of BEAVRS full-core model. (a) Power distribution; (b) Relative standard deviation.

The variance underestimation ratio distribution in Fig. 13 was obtained by comparing the RSDs in a single MC simulation with the reference results. The VURs in each assembly varied from 72% to 91%. The VURs of the assemblies close to the center were lower than the average. By contrast, the outer assemblies had a VUR above 90%. The average VUR was 87.5%.

Fig. 13Full-screen View

(Color online) RSDs and VURs distribution of a single MC simulation in BEAVR full-core model. (a) Relative standard deviation; (b) VURs.

4.2.2

Estimation procedure

A batch size of 50,000 was analyzed to observe the FSD inter-cycle correlation length. Fig. 14 shows that the FSD average SW distance evolves with the interval cycle length. A visual check of the curves showed that the FSD inter-cycle correlation length of this model was approximately 500.

Fig. 14Full-screen View

FSD average SW distance evolves with interval cycle length in BEAVR full-core model

The initial FSD average SW distance was approximately 1/10 of the stationary SW distance. Thus, the variance underestimation problem is more severe in this model than in the sphere array model. This conclusion is consistent with the results presented in Figs. 7 and 13.

4.2.3

Elimination procedure

We set different batch lengths to verify that the FSD inter-cycle correlation length was the optimum batch length. The batch lengths are listed in Table 2. In the SeVUE method, only the case where $l_{\text{b}}=500$ should be calculated.

Table 2Full-screen View

Parameters of batch methods for BEAVRS full-core model

Serial number	Batch length/cycle	Number of batches	Number of active cycles
1	1	20,000	20,000
2	2	10,000	20,000
3	10	2000	20,000
4	20	1000	20,000
5	40	500	20,000
6	80	250	20,000
7	100	200	20,000
8	125	160	20,000
9	200	100	20,000
10	250	80	20,000
11	400	50	20,000
12	500	40	20,000
13	625	32	20,000
14	800	25	20,000
15	1000	20	20,000

Show more

Figure 15 shows that the average VURs evolve with the batch length. The average VUR was approximately 87% when the batch length was one cycle. It gradually declined with the batch length increases. When the batch length exceeded 500, the average VURs were within $\pm 5\%,$ and the curve fluctuated. Thus, 500 is the optimum batch length for this model.

Fig. 15Full-screen View

Average VUR evolves with batch length in BEAVR full-core model

The results of the single MC simulation with $l_{\text{b}}=500$ are shown in Fig. 16. The VUR of each assembly varies from -10.5% to 13.4%, with an average VUR of -1.3%. The differences were acceptable, considering the fluctuations in the reference standard deviations and batch method results.

Fig. 16Full-screen View

(Color online) Relative standard deviation and VURs distributions of an MC simulation with $l_{\text{b}}=500$. (a) Relative standard deviation; (b) VURs

Discussion

Relatively large tallies were used to exemplify the efficiency of the SeVUE method, and this method effectively eliminated the variance underestimation problem in MC criticality simulations. The FSD inter-cycle correlation length was validated to be independent of the batch size. The ratio between the initial FSD average SW distance and stationary SW distance in the Estimation procedure can be a potential reference for the variance underestimation of large local tallies.

The penalty is that an additional Estimation calculation must be conducted to estimate the FSD inter-cycle correlation length of every model. In addition, the FSD inter-cycle correlation length may be hundreds of cycles in the reactor full-core model. Thousands of active cycles and billions of neutron histories are required to eliminate the variance underestimation problem. Moreover, the standard deviation of the tally results may have $\pm 10\text{\%}$ fluctuations in the VUR when a proper batch length is set.