Critical Buckling Generation of TCA Benchmark by the B1 Theory-Augmented Monte Carlo Calculation and Estimation of Uncertainties

The Korea Atomic Energy Research Institute (KAERI) has developed the DeCART2D 2-dimensional (2D) method of characteristics (MOC) transport code and the MASTER nodal diffusion code and has established its own two-step procedure. For design code licensing, KAERI prepared a critical experiment on the verification and validation (V&V) of the DeCART2D code. DeCART2D is able to perform the MOC calculation only for 2D nuclear fuel systems, such as the fuel assembly. Therefore, critical buckling in the vertical direction is essential for comparison between the results of experiments and DeCART2D. In this study, the B1 theory-augmented Monte Carlo (MC) method was adopted for the generation of critical buckling. To examine critical buckling using the B1 theoryaugmented MC method, TCA critical experiment benchmark problems were considered. Based on the TCA benchmark results, it was confirmed that the DeCART2D code with the newly-generated critical buckling predicts the criticality very well. In addition, the critical buckling generated by the B1 theory-augmented MC method was bound to uncertainties. Therefore, utilizing basic equations (e.g., SNU S/U formulation) linking input uncertainties to output uncertainties, a new formulation to estimate the uncertainties of the newly generated critical buckling was derived. This was then used to compute the uncertainties of the critical buckling for a TCA critical experiment, under the assumption that nuclear cross-section data have uncertainties.


Introduction
In nuclear design and analysis, a conventional two-step procedure, involving assemblywise lattice transport and whole core nodal calculations, has been widely and efficiently utilized. This two-step procedure involves unavoidable approximations, such as energy group condensation and assembly-wise homogenization. Energy group condensation and assembly-wise homogenization can be used to establish few group constants (FGCs) from continuous or ultra-fine energy group cross sections and to simplify cell-wise geometric information data. These assembly-wise FGCs are utilized to conduct whole-core nodal calculations. Moreover, a correction of the leakage spectrum (e.g., the B 1 method) is applied to the whole-core calculations for improved accuracy. As a part of the localization of nuclear reactor development in Korea, the Korea Atomic Energy Research Institute (KAERI) developed the DeCART2D [1] 2-dimensional (2D) method of characteristics (MOC) transport code and the MASTER [2] nodal diffusion code and established its own two-step procedure. The DeCART2D code generates assembly-wise FGCs and pin-to-box factors (i.e., pin-wise power distribution), which are used in the MASTER nodal whole-core calculations. To apply the DeCART2D/MASTER code system to nuclear design and analyses, KAERI is now in the progress of licensing the nuclear design analysis code system. Accordingly, KAERI prepared a critical experiment for verification and validation (V&V) of the DeCART2D code. The main purpose of this critical experiment was to check whether the DeCART2D code could accurately provide the pin-to-box factor.
The DeCART2D code can perform the MOC calculation only for a 2D based nuclear system. Therefore, to compare the pin-to-box factors from the critical experiment with those yielded by the DeCART2D code, critical buckling and diffusion coefficients in the vertical direction are necessary. In a critical experiment with a finite cylinder system, vertical buckling generally can be obtained using the extrapolation distance from the activation measurements and the active height of the core.
As an alternative way to obtain the vertical critical buckling, the Monte Carlo (MC) method [3][4][5] is an attractive approach because it provides highly accurate and efficient solutions in the field of nuclear reactor design and analysis. Yamamoto [6,7] developed the heterogeneous B 1 method by introducing a complex-value weight into MC particle transport simulations. The complex-value-weight-based heterogeneous B 1 method can provide directional diffusion coefficients and critical buckling based on leakage-corrected k eff -eigenvalue calculations. Meanwhile, the B 1 theory-augmented MC method can be used to generate fuel pin cell or fuel assembly FA-homogenized few-group diffusion theory constants (FGCs) in the same manner as deterministic-based diffusion theory computations of power reactors [8][9][10][11]. In previous studies [8][9][10], we presented the B 1 theory-augmented MC method for generating homogenized FGCs including diffusion coefficients, critical buckling, and migration area, and verified the results through a nuclear core design analysis of Yonggwang Nuclear Unit 4.
In general, a result yielded by an MC calculation always includes stochastic uncertainty, because the MC method uses a finite number of samples (e.g., neutron particles) through the use of random numbers. Accordingly, the critical buckling by the B 1 theoryaugmented MC method may be bound partly to the statistical uncertainties inherent in the MC method, and partly to the uncertainties of the input data such as the cross section, geometric data, and material composition. In previous studies, we introduced basic mathematical equations aimed at quantifying the uncertainties of FGCs over burnup. These basic mathematical equations from the so-called 'SNU S/U formulation' [12] can be utilized to link input uncertainties to output uncertainties. In the same manner as the uncertainties of FGCs, the SNU S/U formulation can be easily adapted to quantify the uncertainties of the critical buckling generated by the B 1 theory-augmented MC method.
The purpose of this study was to generate critical buckling, based on the B 1 theoryaugmented MC method, for DeCART2D calculations and to quantify the uncertainties of the critical buckling due to stochastic uncertainties and the uncertainties of the input parameters (e.g., cross section, geometric data, and material composition). In this study, the critical buckling for the tank-type critical assembly (TCA) benchmark problem [13,14], which provides the measured values for the critical water level and extrapolation for each case distance, is generated using the B 1 theory-augmented MC method. Furthermore, the uncertainties of the critical buckling are estimated using the SNU S/U formulation. Section 2 briefly describes the generation of the critical buckling using the B 1 theoryaugmented MC method and the direct fitting method. Section 3 explains how to quantify the uncertainties of the critical buckling via the B 1 theory-augmented MC method using the SNU S/U formulation. In Section 4, we apply the newly implemented uncertainty quantification capability to perform the TCA critical facility benchmark for the generation of the critical buckling. In this section, we present the uncertainties in the critical buckling that occur due to the uncertainties in nuclear cross section data and the tolerance of the geometric and material composition data. The conclusions and summations are presented in Section 5. This section briefly explains how to obtain the critical buckling via the B 1 theoryaugmented MC method. The conventional multi-group B 1 equations in a fine group representation are as below: Here, g and g refer to the energy group. Σ t,g and Σ f ,g are the macroscopic total and fission cross-section, respectively. Σ n gg (n = 0,1) indicates the zero and first-order moment of the double differential scattering cross-section. The other notations are standard and φ g is normalized to satisfy k = ∑ νΣ f ,g φ g , where k is the multiplication factor. The macroscopic cross-section is expressed by where ϕ g is the infinite medium spectrum (IMS). N m j and r m α,g,j indicate the number density and α-type microscopic reaction rate for nuclide j in region m. The microscopic reaction rate and IMS are calculated using the microscopic cross-section or angular neutron flux, as given below: where V m is the volume of the region m, whereas ∆E g is the energy interval of group g neutrons. A conventional MC code tallies the microscopic reaction rates and IMS for the total, fission, and scattering macroscopic cross section calculations using the track length and collision estimator method. Using the macroscopic cross sections, the multigroup B 1 equations, Equation (1), are solved through iterative adjustments of B 2 until the corresponding k = 1. The solution φ g and B 2 satisfying k = 1 are the so-called critical spectrum (CS), φ B g , and the critical buckling, B c .

Critical Buckling Generation Using the Fitted Cosine Function
Based on an experiment, the axial buckling [15] for a critical condition can be calculated using the extrapolated length from its activation measurements. In a finite cylinder reactor, normalized flux distributions for the vertical axis can be described well using a sinusoid with the extrapolated length. If a finite cylinder reactor has an extrapolated height H and is centered about z = 0, the vertical flux distribution and geometric buckling can be written as follows: where H is the height of the core in a finite cylinder reactor and λ z is the extrapolated length or distance. In the same manner as the experiment, it is possible to obtain the axial critical buckling from the axial flux distribution, φ MC , which is calculated using the 3-dimensional MC calculation for the critical condition. The vertical extrapolation length can be determined through a fitting procedure for a cosine function. The fitting method is applied to minimize the sum of the squares of the difference between φ MC and the value (φ FIT ) from the fitted cosine function at the same height position. An R-squared (R 2 ) value is used to quantify how close the MC-calculated values are to the fitted regression cosine function. Equation (8) shows the common definition of R 2 between φ MC and φ FIT .
where N and i indicate the total number of mesh tallies and the mesh index number, respectively. Figure 1 summarizes the procedure used to obtain the critical buckling in the B 1 theory-augmented MC method. The whole procedure is divided into three steps. The uncertainty of the critical buckling is calculated step-by-step, because the outputs of any given step become inputs of the next step. The basic equations linking the uncertainties of the inputs to those of the outputs were derived in a previous study [12]. Suppose that Q represents the output parameters at any given step and u i (i = 1, 2, · · · , I) is a set of uncertain input data for Q calculation. Formally, Q may be written as follows: cov [ , ] ii uu is the covariance of two uncertain input variables, ui' and ui''. The partial de-147 rivatives in Equation (12) are calculated using the following approximation with the direct 148 subtraction method or the MC perturbation technique.

Estimation of the Critical Buckling Uncertainties Using the SNU S/U Formulation
150 Figure 1. Flowchart of critical buckling generation using the B1 theory-augmented MC method. 151 For the application of Equations (11) and (12) to uncertainty quantifications of any 152 one-step calculations in critical buckling generation, a functional relation between the un-153 certain input/output (I/O) variables such as ui and Q must be established. Table 1 shows 154 the functional relation between I/O variables for each step in the critical buckling genera-155 tion using the B1 theory-augmented MC method. According to the SNU S/U formulation, 156

STEP III Fine Group B1 Calculations Number Density with Uncertainty
Cross Section and its Covariance Data

Fine Group Constants and its Uncertainty
Critical Buckling and its uncertainty Through the deviation of the SNU S/U formulation, linking uncertainties using the 1st order Taylor expansions [8], the uncertainty of Q is expressed by where the angular bracket indicates the expected value of quantity within it and cov[u i , u i ] is the covariance of two uncertain input variables, u i and u i" . The partial derivatives in Equation (12) are calculated using the following approximation with the direct subtraction method or the MC perturbation technique.
For the application of Equations (11) and (12) to uncertainty quantifications of any onestep calculations in critical buckling generation, a functional relation between the uncertain input/output (I/O) variables such as u i and Q must be established. Table 1 shows the functional relation between I/O variables for each step in the critical buckling generation using the B 1 theory-augmented MC method. According to the SNU S/U formulation, the variance of the critical buckling from the B 1 calculations at step III, σ 2 (B c ), is formulated by using Equation (10) as follows: Table 1. Functional relation between input and output variables for each step in critical buckling generation using the B 1 theory-augmented MC method. Step

Functional Relation Between the I/O Variables
Energies 2021, 14, 2578 Because there is no statistical uncertainty in solving B 1 multi-group equations, σ 2 STATS (B c ) = 0, using the functional relation between I/O variables in step II, one can obtain the covariance term between macroscopic cross sections, cov[Σ α ,g , Σ α ,g ], in terms of new covariance terms involving the uncertain input parameters, namely, number density, microscopic reaction rate, and IMS, as shown below: In the same manner, the unknown covariance terms in Equation (15) can be derived using the functional relations between I/O variables in step I. Equation (16) through (20) indicate these covariance terms. They can be expressed in terms of covariance terms involving the uncertain input parameters, such as microscopic cross section and number density. The partial derivatives are calculated using the direct subtraction method, as shown in Equation (13).
Energies 2021, 14, 2578 Note that an examination of the uncertainty of the critical buckling, as shown in Equation (14), ultimately entails calculations of the covariance between number densities or microscopic cross sections. In this study, only the covariance between cross sections is considered because the covariance terms directly related to number density are set at 0. The covariance between microscopic cross sections, cov[σ α ,g ,i , σ α ,g ,i ], can be obtained by processing the covariance files of the evaluated nuclear data library files using the ERRORR module of NJOY [16].

Estimation of the Critical Buckling Uncertainties using the Stochastic Sampling Method
The stochastic sampling (S.S.) method, which is often called the brute force method, provides another approach to determining the uncertainties of the critical buckling. The mean of the uncertain input parameter, u i , and the covariance between u i and u j uncertain input data are defined by Supposing that C u is the covariance matrix defined by cov[u i , u j ] and that a lower triangular matrix B is known through the Cholesky decomposition of C u , we then have where B T is the transpose matrix of B. Then, if C u is symmetrical and positive definite, one can obtain a sample set by where X is the mean vector defined by the mean values from Equation (21) and Z is a random normal vector calculated directly from a random sampling of the standard normal distribution using the Box-Muller method. X κ indicates the sample set generated from a κ-th sampled Z vector.

TCA Benchmark Problem
The TCA benchmarks [13,14] are critical experiment series for light-water moderated lattices using 17 × 17 2.6 w/o UO 2 fuel rods arranged in a square array. The fuel rod array is surrounded by light-water in a 1.8 m core tank, whereas the pin-pitch and fuel pellet radii are 1.956 cm and 0.625 cm. Basically, the critical sizes were determined by adjusting the light water level. The TCA benchmark provides various measured parameters, such as the temperature coefficient, the extrapolation length, and the critical buckling. The experiments were performed using various configurations between 1963 and 1975. Based on these experiments, eighteen configurations were provided in the International Criticality Safety Benchmark Problem (ICSBEP) [13]. Among them, five experimental cases (i.e., case 4-case 8) were selected for this study. Table 2 shows the specifications of the five selected benchmark problems. Figure 2a extrapolation distance (λ z ) are obtained from the experimental results. Using the measured extrapolation distances, the total critical buckling can be calculated as where P A is the pitch of a fuel rod array.

TCA Benchmark Problem
The TCA benchmarks [13,14] are critical experiment series for light-water moderated lattices using 17 × 17 2.6 w/o UO2 fuel rods arranged in a square array. The fuel rod array is surrounded by light-water in a 1.8 m core tank, whereas the pin-pitch and fuel pellet radii are 1.956 cm and 0.625 cm. Basically, the critical sizes were determined by adjusting the light water level. The TCA benchmark provides various measured parameters, such as the temperature coefficient, the extrapolation length, and the critical buckling. The experiments were performed using various configurations between 1963 and 1975. Based on these experiments, eighteen configurations were provided in the International Criticality Safety Benchmark Problem (ICSBEP) [13]. Among them, five experimental cases (i.e., case 4-case 8) were selected for this study. Table 2 shows the specifications of the five selected benchmark problems. Figure 2a,b display the horizontal and vertical cross sections for the TCA benchmark problems. The 1.8 m water core tank of the critical facility was simplified as a 30-cm-thick side water reflector in the benchmark model. For each problem, the critical water level (Hc), measured horizontal extrapolation distance ( h  ), and vertical extrapolation distance ( z  ) are obtained from the experimental results. Using the measured extrapolation distances, the total critical buckling can be calculated as where PA is the pitch of a fuel rod array.

Horizontal Plane (X-Y)
Number

Critical Buckling Generation for TCA Benchmark Problem
For each TCA benchmark problem, the critical buckling was generated using the McCARD [3] code, via both the B 1 theory-augmented MC method and the direct fitting method. In the B 1 method, the macroscopic cross sections for solving B 1 multi-group equations were generated using MC calculations based on 50 inactive and 100 active cycles with 100,000 histories per cycle. The ACE-format continuous energy cross section libraries were processed based on the ENDF/B-VII.1 evaluated nuclear data library. The 47-group structure [1] was adopted for the B 1 multi-group calculations. In the direct fitting method, the vertical extrapolation lengths were calculated using the MC calculations and the horizontal extrapolation lengths were taken from [13]. For the calculation of vertical extrapolation length, a single fuel rod located at the center of a fuel array was axially divided into one hundred meshes, as shown in Figure 3. Figure 4 shows the normalized vertical flux distributions and their fitted cosine functions for each problem. Table 3 shows the vertical extrapolation lengths determined using the direct fitting method. For all problems, the R 2 values for their extrapolation lengths were more than 0.999. The maximum difference in the vertical extrapolation length between the experiment and the calculation was around 5.0% for problem 5. were processed based on the ENDF/B-VII.1 evaluated nuclear data library. The 47-group structure [1] was adopted for the B1 multi-group calculations. In the direct fitting method, the vertical extrapolation lengths were calculated using the MC calculations and the horizontal extrapolation lengths were taken from [13]. For the calculation of vertical extrapolation length, a single fuel rod located at the center of a fuel array was axially divided into one hundred meshes, as shown in Figure 3. Figure 4 shows the normalized vertical flux distributions and their fitted cosine functions for each problem. Table 3 shows the vertical extrapolation lengths determined using the direct fitting method. For all problems, the R 2 values for their extrapolation lengths were more than 0.999. The maximum difference in the vertical extrapolation length between the experiment and the calculation was around 5.0% for problem 5.     Table 4 shows the critical buckling of McCARD, calculated using the B1 method and 238 the direct fitting method. The relative differences in critical buckling between the refer-239 ence and McCARD range from 0.1% to 1.1%. The maximum relative difference in the crit-240    Table 4 shows the critical buckling of McCARD, calculated using the B 1 method and the direct fitting method. The relative differences in critical buckling between the reference and McCARD range from 0.1% to 1.1%. The maximum relative difference in the critical buckling occurred for problem 1. In all problems, the B 1 calculations were conducted under the equivalent geometric data because the B 1 method uses an infinite medium bounded by surfaces with a reflective boundary condition. Accordingly, the critical buckling results determined by the B 1 method were the same, as shown in Table 4. To confirm the accuracy of the criticality prediction, DeCART2D lattice code calculations were performed using the critical buckling. For the DeCART2D calculation, a ray spacing of 0.02 cm and two polar angles of 90 • were used as the ray tracing option. All the DeCART2D calculations were conducted with the ENDF/B-VII.1 based 47-group cross section library [17]. Table 5 shows the effective multiplication factors (k eff ) of DeCART2D for the five TCA benchmark problems. The maximum differences in k eff for the reference is about 143 pcm whereas that for the B 1 method is about 187 pcm. The direct fitting method gives a maximum error of 217 pcm. It is well known that the general design review criterion (DRC) of critical boron concentration (CBC) in a typical PWR start-up and operation is 500 pcm. Therefore, it can be confirmed that all DeCART2D results agree with the experiments within the DRC of CBC. Table 5. k eff of DeCART2D with critical buckling from vertical extrapolation length.

Problem
No.

Uncertainty Analysis of Critical Buckling Based on Cross Section Uncertainties
The SNU S/U formulations to quantify the uncertainties arising from the statistical uncertainty and the uncertain input data, which were derived in Section 3.1, were incorporated into the B 1 method-based critical buckling generation module of the McCARD code. The module was then used to estimate the uncertainties of the critical buckling for the TCA benchmark problems. The nuclear cross-section data, including the covariance data files for computing both the critical buckling and their uncertainties, were obtained from the ENDF/B-VII.1 and ENDF/B-VIII.0 evaluated nuclear data libraries. The covariance data files of only the two major uranium isotopes (i.e., 235 U and 238 U) and 16 O were considered. Table 6 presents the percentile relative error (%RE), which is defined as 100 × σ(B 2 c )/B 2 c , for each critical buckling. All relative errors of critical buckling caused only by statistical uncertainties were less than 0.11% whereas those from both statistical and microscopic cross section uncertainties for ENDF/B-VII.1 and ENDF/B-VIII.0 were less than 3.18% and 1.49%, respectively. It is clear that the most significant contributor of both evaluated nuclear data file libraries to σ(B 2 c ) is the uncertainties of the 238 U capture cross section (MT = 102).

291
The S.S. method can be easily utilized to quantify uncertainties of nuclear design pa-292 rameters, including critical buckling, through calculations of the inputs, which are sam-293 pled from the standard deviation or tolerance of the input parameters. The TCA experi-294 ment benchmark provides the uncertainty of input design data such as the uranium com-295 position, the diameter of the UO2 pellet, the inner and outer diameter of the aluminum 296 alloy cladding, the pin pitch, and the density of the UO2 pellet. In this section, the critical 297 buckling uncertainty based on these uncertainties was estimated using the S.S. method. 298 All S.S. calculations were performed using McCARD and MIG (McCARD Input Genera-299 tor) code [9]. For each case, the 100 McCARD calculations were conducted using the 300 ENDF/B-VII.1 cross section library. Table 7 presents the uncertainties of the critical buck-301

Uncertainty Analysis of Critical Buckling Based on Input Design Data
The S.S. method can be easily utilized to quantify uncertainties of nuclear design parameters, including critical buckling, through calculations of the inputs, which are sampled from the standard deviation or tolerance of the input parameters. The TCA experiment benchmark provides the uncertainty of input design data such as the uranium composition, the diameter of the UO 2 pellet, the inner and outer diameter of the aluminum alloy cladding, the pin pitch, and the density of the UO 2 pellet. In this section, the critical buckling uncertainty based on these uncertainties was estimated using the S.S. method. All S.S. calculations were performed using McCARD and MIG (McCARD Input Generator) code [9]. For each case, the 100 McCARD calculations were conducted using the ENDF/B-VII.1 cross section library. Table 7 presents the uncertainties of the critical buckling from geometric and material data uncertainties for Problem 1. The percentile relative errors (%RE) of the critical buckling by the UO 2 pellet diameter uncertainty, the aluminum cladding inner and outer diameter uncertainty, and the pin pitch uncertainty are 0.41%, 0.04%, 0.07%, and 0.12%, respectively. The uncertainties of the critical buckling by the 235 U enrichment uncertainty and the UO 2 pellet density uncertainty are respectively 0.15% and 0.28%. Figure 6 displays the distributions of critical buckling from each input design data sampling. It was observed that the critical buckling increases as the fuel pellet diameter, pin pitch, fuel pellet number density, or 235 U enrichment increases, whereas it slightly decreases as the cladding inner or outer diameter increases. These linear relationships for each case depend on the change in the reactivity. The increase in the fuel pellet diameter leads to an increase in the reactivity, the leakage term for the critical state, and the critical buckling, in sequence.

Conclusions
In this study, the B 1 theory-augmented Monte Carlo (MC) method was adopted for the generation of critical buckling. To examine the critical buckling using the B 1 theoryaugmented MC method, five TCA critical experiment benchmark cases were considered. For each case, the effective multiplication factor of DeCART2D, which was calculated with the critical buckling determined by the B 1 method, agreed with the experimental value within the well-known general criteria of CBC in a typical PWR start-up and operation. Therefore, it can be concluded that the B 1 theory-augmented MC method can provide accurate and reliable critical buckling for a practical system.
Moreover, a new formulation to estimate the uncertainties of the critical buckling as determined by the B1 theory-augmented MC method was derived from the so-called 'SNU S/U formulation', and this can be used to link input uncertainties to output uncertainties. The SNU S/U formulation was able to explain the uncertainty propagation in the critical buckling generation. It was then utilized to compute the uncertainties of the critical buckling for a TCA critical experiment under the assumption that nuclear cross-section data have uncertainties. This application of the 'SNU S/U formulation' also shows that this formulation can serve as an effective and useful tool to examine the influence of nuclear data uncertainty on nuclear core design and safety. Moreover, it can be utilized to identify the significant contributors to the uncertainties of the other nuclear core design parameters.
The procedure involving the new formulation will be particularly useful as a practical way to provide the uncertainties of the critical buckling required for V&V of the DeCART2D code.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.