PSI Methodologies for Nuclear Data Uncertainty Propagation with CASMO-5M and MCNPX: Results for OECD/NEA UAM Benchmark Phase I

Abstract

Capabilities for uncertainty quantification (UQ) with respect to nuclear data have been developed at PSI in the recent years and applied to the UAM benchmark. The guiding principle for the PSI UQ development has been to implement nonintrusive “black box” UQ techniques in state-of-the-art, production-quality codes used already for routine analyses. Two complimentary UQ techniques have been developed thus far: (i) direct perturbation (DP) and (ii) stochastic sampling (SS). The DP technique is, first and foremost, a robust and versatile sensitivity coefficient calculation, applicable to all types of input and output. Using standard uncertainty propagation, the sensitivity coefficients are folded with variance/covariance matrices (VCMs) leading to a local first-order UQ method. The complementary SS technique samples uncertain inputs according to their joint probability distributions and provides a global, all-order UQ method. This paper describes both DP and SS implemented in the lattice physics code CASMO-5MX (a special PSI-modified version of CASMO-5M) and a preliminary SS technique implemented in MCNPX, routinely used in criticality safety and fluence analyses. Results are presented for the UAM benchmark exercises I-1 (cell) and I-2 (assembly).

1. Introduction

The OECD/NEA benchmark for uncertainty analysis in modeling (UAM) was launched a few years ago to promote the development, assessment, and integration of comprehensive uncertainty quantification (UQ) methods in best-estimate multiphysics coupled simulations of LWRs during normal as well as transient conditions [1]. Although very ambitious by nature (due to the complexity of the task to treat all potential sources of uncertainties), the benchmark has nevertheless achieved one of its first objectives, namely, to constitute a major (if not the main) international framework to drive forward the development of methodologies for the propagation of nuclear data uncertainties in reactor simulations. This topic was proposed as the first phase of the benchmark, and since research in precisely this area was at the same time being launched within the STARS project [2] at the Paul Scherrer Institut (PSI), participation to this benchmark was considered as a timely and highly valuable opportunity to complement the development and assessment of the PSI methods. In that context, two parallel lines of development were in fact initiated at PSI. On the one hand, the development of a UQ methodology for the propagation of neutronic uncertainties in the deterministic CASMO/SIMULATE/SIMULATE-3 K chain of reactor analysis codes and used for safety assessment of the Swiss reactors was launched. On the other hand, the development of a corresponding UQ methodology for neutron transport simulations with the stochastic continuous-energy MCNPX and with primary emphasis on criticality safety was recently initiated. In this paper, the principles and concepts of both methodologies are first summarized. Then, the results obtained for the UAM Phase I benchmark cases analyzed so far are presented. The primary focus is given to the CASMO-5 M analyses conducted so far for Phase I-1, aimed at cell physics, and Phase I-2, dedicated to lattice physics. Regarding MCNPX, the first set of solutions obtained for the hot-zero-power pin cell cases of Phase I-1 will also be presented.

1.1. Motivation

In order to rigorously establish the accuracy (or bias) of the so-called best-estimate codes, the precision (or uncertainty) must be quantified. (The measure of accuracy is bias: low accuracy implies a large bias and high accuracy implies a small bias. The measure of precision is uncertainty: low precision implies large uncertainty and high precision implies small uncertainty.) This includes propagation of input uncertainty (all inputs are really distributions) to output uncertainty, which is the basic task of UQ. The most straightforward benefit of UQ is the new information about the distribution of outputs which can be used to qualify designs and/or provide confidence in results. However, with UQ a much more rigorous validation procedure is also available and the value of this should not be underestimated. With UQ, one can compare calculations with uncertainty 𝐶±𝜎_𝐶 to experimental results with uncertainty 𝐸±𝜎_𝐸 using overlap testing, instead of the conservative assumption of no uncertainty in calculations (𝜎_𝐶=0) or subjective use of expert judgment to decide if 𝐶 is close enough to 𝐸±𝜎_𝐸. A best-estimate code (and its validation) should avoid such conservative assumptions and expert judgment by definition.

1.2. Preliminaries

Consider input, 𝑥, and output, 𝑦, with nominal values 𝑥₀ and 𝑦₀ and perturbed values 𝑥^′ and 𝑦^′. When sampling perturbed values 𝑥^′ from a distribution, the 𝑛th sample of the input is 𝑥^′(𝑛) and the corresponding output is 𝑦^′(𝑛). With computer codes, there is typically a large set of input and output, which may be denoted 𝐱₀ and 𝐲_𝟎 for nominal sets and 𝐱^′, and 𝐲^′ for perturbed sets.

1.3. First-Order UQ Using Uncertainty Propagation

The cornerstone of local, first-order UQ methods is the capability to calculate sensitivity coefficients: 𝑆≡𝜕𝑦𝜕𝑥_𝑥=𝑥0𝑥₀𝑦₀, (1) which are vital to sensitivity analysis (SA). It is very convenient to introduce the definition of a perturbation factor, 𝑝, such that the perturbed input 𝑥^′=𝑝𝑥₀ and the corresponding output factor, 𝑞, 𝑦^′=𝑞𝑦₀. Thus the sensitivity coefficient may be written simply as 𝑆≡𝜕𝑞𝜕𝑝_𝑝=1; (2)

Nonintrusive SA can then be implemented simply as a numerical differentiation of 𝑞 with respect to p, referred to here as direct perturbation (DP) as in [3]. Two factors make it difficult to use DP in an automated manner to obtain accurate estimates of 𝑆: (i) due to finite-precision arithmetic, “too small” perturbations do not change the output significantly and (ii) due to unknown relationships between inputs and outputs, “too small” perturbations for one input may be “too large” for another. With low-order numerical differentiation formulas (e.g., first-order finite differences) especially, “too large” perturbation can greatly increase the approximation error. Due to the relatively high cost of calculations in nuclear simulations, low-order formulas are typically used.

Although DP can be straightforwardly extended to simultaneously estimate 𝑆 for multiple outputs, it cannot handle simultaneous input perturbations; that is, input perturbations are always one at a time. Thus with many more input parameters than output parameters, DP is not very efficient. For nuclear data uncertainty propagation, this was the basic reason behind the development of very efficient (but intrusive) perturbation theory-based algorithms for sensitivity coefficient estimation, for example, in the SCALE code system [4].

Using the calculated sensitivity coefficients in UQ simply requires the classic first-order uncertainty propagation formula [5], shown below for multiple input and output parameters: 𝐕_𝐘=𝐒^𝐓𝐕_𝐗𝐒, (3) with relative variance/covariance matrix (VCM) of the outputs 𝐕_𝐘 in terms of the relative VCM of the inputs 𝐕_𝐗 and the sensitivity coefficients, 𝐒, now a matrix defined as 𝑆_𝑖𝑗=𝜕𝑞_𝑗/𝜕𝑝_𝑖 for input parameter index (row) 𝑖 and output parameter index (column) 𝑗.

1.4. Sampling-Based UQ

Sampling-based UQ, or stochastic sampling (SS), has been historically used for nonlinear systems with few correlated parameters [6]. However, currently SS is increasingly applied for all types of UQ, including neutronics, due to its nonintrusive nature [7], flexibility to handle many uncertain parameters [8], theory for non-parametric tolerance intervals (i.e., Wilks’ formula), and global sampling of the solution space. In order to implement an SS method, one simply has to sample inputs from their distributions (choosing appropriate distributions is another matter), run the code with each sampled set, and analyze the distribution of outputs.

1.4.1. Simple Random Sampling from Multivariate Gaussian Distributions

In the case of the distribution of nuclear data, one generally assumes that the input 𝐱 obeys an 𝑀-dimensional Gaussian (normal) distribution [5]: 𝑓_𝑋(𝐱)=1(2𝜋)^𝑀/2𝐕_𝐗^1/2exp−𝐱−𝐱_𝟎𝐕_𝐗^−𝟏𝐱−𝐱_𝟎^𝑇2, (4) with input VCM 𝐕_𝐗 of dimension 𝑀×𝑀, |𝐕_𝐗| the determinant, and mean 𝐱_𝟎 of dimension 1×𝑀. A matrix of 𝑁 simple random samples 𝐗^′ (𝑁×𝑀) which respects the correlations of data may be constructed as described in [9].(1)Decompose VCM 𝐕_𝐗 using a “Cholesky-like” decomposition (see below), 𝐓^𝐓𝐓=𝐕, where 𝐓 is 𝑀^∗×𝑀.(2)Make 𝑁×𝑀^∗ random samples from the standard normal distribution (zero mean and unit variance) and store the results in the 𝑁×𝑀^∗ matrix 𝐙. (3)The random samples are then given as 𝐱^′(𝐧)=𝐱_𝟎+𝐙^(𝐧)𝐓, where 𝐙^(𝐧) is the 𝑛th row of 𝐙.

The term “Cholesky-like” is used because a true Cholesky factorization requires a (square) symmetric positive definite (SPD) matrix whereas a general VCM can be symmetric positive semidefinite (SPSD), for example, due to perfect (anti) correlation of parameters. In this case the matrix is rank-deficient with rank 𝑀^∗ and 𝐓 is rectangular, with 𝑀^∗<𝑀, and (4) must use the generalized inverse and pseudodeterminant.

1.4.2. Nonparametric Statistics and Wilks’ Formula

Given 𝑁 random samples of a quantity, the formula for the tolerance interval in terms of coverage (a) and confidence (b), without assuming a particular distribution, is known colloquially as Wilks’ formula, due to the seminal work of S. S. Wilks in nonparametric statistics [10]. Nonparametric (or order) statistics is the name given to the set of statistical techniques which does not require data belonging to a particular distribution (e.g. normal) and frequently requires ordering samples, for instance, from least to greatest. For a more complete discussion of nonparametric statistics applied to neutronics calculations, see [8]. In order for Wilks’ formula to be valid, a simple random sampling process must be used; that is, stratified sampling or variance reduction is not allowed according to the theory. For example, with 𝑁=93 samples, a two-sided tolerance interval can be declared as [𝑦_min,𝑦_max], where 𝑦_min and 𝑦_max are the minimum and maximum results from the 93 samples. Such a tolerance interval is guaranteed to contain the (middle) 𝑎=95% of the distribution with 𝑏=95% confidence. Note that, with 𝑎=𝑏, the behavior of log (𝑁) versus log (1−𝑎) is roughly linear.

1.4.3. Sample Statistics

In neutronics UQ, the variance (or standard deviation) is used most often as the measure of uncertainty. With UQ methods based on the uncertainty propagation formula (e.g., DP), the variance of outputs is simply the diagonal of the output VCM. With SS, it is convenient to use the sample variance from sample statistics: 1𝑣=𝑁^𝑁_𝑛=1𝑦^′(𝑛)−𝑦₀², (5) where 𝑦^′(𝑛) is the (perturbed) result of sample 𝑛, 𝑦₀ is the nominal calculation value, and 𝑁 is the total number of samples. It is well known in statistics that the sample variance of a normal distribution is a scaled chi-square distribution of 𝑁−1 degrees of freedom which can be used to provide 𝑁-dependent bounds on the sample variances as shown in [6].

2. Methodology

Although both direct perturbation (DP) and stochastic sampling (SS) schemes are “nonintrusive” by nature, in order to develop UQ techniques for the CASMO-5 M lattice physics code, some source modifications were necessary as CASMO-5M’s nuclear data library is stored in a proprietary binary format and “perturbed libraries” could not be easily created.

For the relatively newer developments concerning UQ with MCNPX, ACE format libraries may be created directly and thus no source code modifications of MCNPX are required. The following sections will first describe the CASMO-5MX code, then the DP and SS techniques as designed for use with CASMO-5MX, and finally the SS technique development for MCNPX.

2.1. General Development of CASMO-5MX

The capability to perturb the nuclear data library of the lattice physics code is the first step in order to perform any “nonintrusive” UQ with respect to nuclear data. Because of the aforementioned proprietary nature of CASMO-5M’s 586-group ENDF/B-VII.R0-based nuclear data library, source code modifications were the easiest way to gain access to this library to perform perturbations. For this purpose, a special module called “PERTXS” and a corresponding cross-section (XS) “perturbation file” was developed. The perturbation file can simply be thought of as a new (optional) input file that demands nuclear data perturbations to apply to the nominal library at runtime. This PSI-modified version of CASMO-5 M will hereon be referred to as CASMO-5MX and the DP technique has been described in [3] and SS technique in [9]. Here, for the reader’s convenience, all necessary elements will be reviewed.

2.1.1. Allowed Nuclear Data Perturbations

Currently CASMO-5MX allows nuclear data perturbations to the following microscopic data for all nuclides in the library (ENDF MT numbers in parentheses).(1)elastic scattering (MT=2),(2)inelastic scattering (MT=4),(3)(𝑛,2𝑛)(MT=16),(4)fission (MT=18), (5)capture (MT=101),(6)average neutrons per fission (MT=452), and(7)average fission spectrum (MT=1018).

In addition, external utilities have been created to perturb any parameter contained in the standard input file, facilitating sensitivity/uncertainty analysis with respect to such parameters as clad thickness, fuel enrichment, and so forth. With nuclear data perturbations, it is important to understand that perturbations are made relative to the existing data on the library; that is, values on the library are not replaced with new perturbed values but increased/decreased by a perturbation factor 𝑝.

2.1.2. Perturbation Formulas

A very convenient feature of CASMO-5MX is that perturbations may be supplied in any group structure, for example, the 19-group “coarse” structure used by default in CASMO-5 M for UO2 assembly method of characteristics (MOC) transport calculations, an arbitrary two-group structure, or the full 586-group library structure. However, using coarse groups for perturbations keeps data files smaller and in most cases, it has been found that using a very fine group structure (e.g., 586 library groups) does not significantly alter the final output uncertainty estimates. (A small study of this will be provided later.) Additionally, because the underlying VCM data is in the SCALE6 44-group structure, it does not make too much sense to go beyond this. Inside CASMO-5MX, the following perturbation formulas are used to map perturbations from the input group structure to the 586-group library structure: 𝑝^𝑔_𝜒=∑_ℎ𝑝^ℎ_𝐼𝜒^ℎ_𝐼𝑤₀^𝑔,ℎ∑_ℎ𝜒^ℎ_𝐼𝑤₀^𝑔,ℎ,𝑝^𝑔_𝜎=∑_ℎ𝑝^𝐻_𝐼𝜎^ℎ_𝐼𝑤₁^𝑔,ℎ∑_ℎ𝜎^ℎ_𝐼𝑤₁^𝑔,ℎ. (6)(7) Equation (6) defines the perturbation factor to be applied to the library fission spectrum in library group 𝑔 with input fission spectrum 𝜒^ℎ_𝐼 and input perturbation 𝑝^ℎ_𝐼, where the input (coarse) energy groups use index ℎ. Equation (7) defines the perturbation factor for a cross-section with standard flux weighting. The weights 𝑤₀^𝑔,ℎ and 𝑤₁^𝑔,ℎ are given by the following equations: 𝑤₀^𝑔,ℎ=∫^{𝐸upper(𝑔,ℎ)}_{𝐸𝑙ower(𝑔,ℎ)}𝑑𝐸∫^𝐸ℎ−1_𝐸ℎ=𝑑𝐸Δ𝐸^𝑔,ℎΔ𝐸^ℎ,𝑤₁^𝑔,ℎ=∫^{𝐸upper(𝑔,ℎ)}_{𝐸lower(𝑔,ℎ)}𝜙(𝐸)𝑑𝐸∫^𝐸𝑔−1_𝐸𝑔𝜙(𝐸).𝑑𝐸 (8)(9) The upper and lower bounds of the numerator integrals are basically the union grid boundaries for the union of group 𝑔 and group ℎ; therefore, the weights are only nonzero where groups overlap. If the supplied perturbation group structure and the library group structure are aligned and there is only one ℎ group for one or more 𝑔 groups (i.e., the input structure is coarser) and the formulas reduce considerably to 𝑝^𝑔=𝑝^ℎ_𝐼 (Figure 1). This means one may perturb the library using only relative information, that is, a set of 𝑝 values. However, when the perturbation group structure is nonaligned with the library, then both the weights and cross-section factors, for example, 𝜎^ℎ_𝐼 and 𝑤₁^𝑔,ℎ, do not cancel in (6) and (7) and dependence on the intragroup weighting functions and cross-sections are introduced. This means that one cannot simply use the perturbations 𝑝 and an approximate spectrum, for example, 𝜙(𝐸)∝1/𝐸, and reference values for the data in the perturbation group structure, for example, 𝜎^ℎ_𝐼, must be provided as well. In CASMO-5MX, the weighting is assumed 1/𝐸 and the reference values in the SCALE6 VCM library are used.

Figure 1 

Aligned and nonaligned coarse perturbation group structures.

The ability to supply perturbations in any group structure effectively gives the user the ability to generate sensitivity profiles at different resolutions for different reactions. For example, to simply evaluate the order-of-magnitude effect for a particular reaction, two-group perturbations could be supplied. If the sensitivity is high, perturbations in a finer structure could be made to generate a refined sensitivity profile. The limiting resolution is simply that of the underlying 586-group library.

2.1.3. Nuclear Data Variance/Covariance Matrices

The uncertainty in groupwise nuclear data is typically expressed only in terms of variance/covariance matrices (VCMs), which implies an underlying Gaussian (normal) distribution of the data. At the single-nuclide, single-reaction level with 𝐺 energy groups, this is a matrix of size 𝐺×𝐺 with the diagonal elements giving the groupwise variance and off-diagonal elements giving the covariance between two groups. Close groups tend to be highly correlated, for example, it is improbable that the data in one fast group would increase and the next one down would decrease. Component cross-sections (e.g., scattering and capture) tend to be anticorrelated, as they must sum to the total cross section. Because measurements are frequently made on compounds, not single nuclides, there is additional correlation between some of the single-nuclide data.

With correlations that cannot be neglected and huge datasets (e.g., 300 nuclides with 44 energy groups and 6 reactions is about 80,000 “inputs”), nuclear data uncertainty propagation is difficult and unique. Because this data is only recently being fully utilized, there are few choices for robust and reasonable VCM evaluations. The SCALE6 VCM [4] data is among the most widely used and developed for these purposes and has been used exclusively in this work, with one single additional approximation due to current limitations in some of the processing tools: cross-nuclide covariances (e.g., Pu²³⁹ fission anticorrelated with Pu²⁴¹ fission) are neglected. The data available on the VCM library and the data which may be perturbed with CASMO-5MX are for the most part consistent. Two exceptions are that the VCM library contains data for each partial capture reaction MT=102–109 and the CASMO-5 M data library combines elastic, inelastic, (𝑛,2𝑛), and (𝑛,3𝑛) into a single “scattering matrix.” The first issue is easily circumvented using the uncertainty propagation formula in (3) to combine the partial VCMs for MT=102 to 109 into a single MT=101 VCM [11]. The second issue dealing with combined scattering is described in the next section.

Because the SCALE6 VCM library is provided in a 44-group structure, nonaligned with the CASMO-5MX library structure, there are two options to use this data:(1)make perturbations in the 44-group structure, relying on (6) and (7) to map these perturbations to 586-groups,(2)convert the 44-group VCMs to a different group structure, ideally to a coarse group structure aligned with the library.

The second option has been investigated and the code ANGELO which performs the conversion has been provided for the purposes of this benchmark [12]. Although its applicability has not been rigorously determined, for converting 44-group SCALE6 VCMs to coarse 8-, 19-, and 31-group structures of CASMO-5MX, the scheme seems reasonable.

2.1.4. Scattering Matrix Perturbations and External Scattering Fraction Data

Many lattice physics codes, including CASMO-5 M, store a single “combined scattering matrix” for each nuclide, lumping elastic and inelastic scattering with the (𝑛,2𝑛) and (𝑛,3𝑛) reactions. Additionally, on the VCM library, uncertainty information for these reactions is only present in “1D” or “vector” form; that is, it has been “summed” over all final energy groups. With these two constraints, perturbations could originally [3] only be applied in the following manner to the combined scattering matrix: 𝜎_𝑠^{′𝑔→𝑔′}=𝑝^𝑔_𝑆𝜎_𝑠^{𝑔→𝑔′}, (10) where perturbation 𝑝^𝑔_𝑆 depends only on the initial group 𝑔 and is applied identically to all final groups, 𝑔^′, in combined scattering matrix, 𝜎_𝑠^{𝑔→𝑔′}. One upside to this type of perturbation is that the mapping formula from (9) can still be used for scattering perturbations. To denote this type of special perturbation, the special MT number MT=13 was introduced to denote “combined scattering” perturbations within the CASMO-5MX system.

However, it became apparent that the combined treatment tends to underestimate the uncertainty due, in particular, to inelastic scattering in U-238 [11], which is actually one of the dominant sources of uncertainty for many responses. However, an approach to separate these effects was described in [9], where one can perform additional NJOY calculations to estimate the so-called “scattering fractions,” that is, fractions of the combined scattering matrix which are due to elastic, inelastic, and so forth. The scattering fractions become an auxiliary library to be used when separation of effects is important. In this case, the scattering matrix perturbation formula becomes 𝜎^{′𝑔→𝑔′}_𝑠=𝑝^𝑔_{𝑆,𝑒𝑙}𝑓^{𝑔→𝑔′}_𝑒𝑙+𝑝^𝑔_{𝑆,𝑖𝑛}𝑓^{𝑔→𝑔′}_𝑖𝑛+𝑝^𝑔_𝑛,2𝑛𝑓^{𝑔→𝑔′}_𝑛,2𝑛+𝑝^𝑔_𝑛,3𝑛𝑓^{𝑔→𝑔′}_𝑛,3𝑛𝜎^{𝑔→𝑔′}_𝑠, (11) where the 𝑓 terms are the scattering fractions, tabulated for each nonzero 𝑔,𝑔^′ pair for that reaction. Currently the scattering fractions have been prepared for U-235 and U-238 only, and only at a temperature of 500 K and a background cross section of 40 barns, after some initial studies, found them to be remarkably constant with respect to temperature and background cross-section variations.

2.1.5. Resonance Self-Shielding

The way that resonance self-shielding is performed in CASMO-5M makes it difficult to perturb nuclear data before the resonance self-shielding calculation. Therefore, the resonance self-shielded and infinitely dilute data are perturbed by the same factor 𝑝, which neglects the effect changes in the data have on self-shielding. Because self-shielding is a “negative” type of feedback, the current approach in CASMO-5MX is thought to produce slightly higher uncertainties, but comparisons to SCALE6 TSUNAMI, which does include the effect, have not shown a significant effect [3, 9]. The difference should be most noticeable with strong and highly uncertain resonances, for which perhaps the U-238 dominated systems tested so far do not qualify.

2.2. Direct Perturbation with CASMO-5MX/DP

The main difficulties applying the DP technique to calculate sensitivity coefficients, namely, fixed-precision and eliminating second-order and higher effects, have been overcome using an adaptive technique [3] in which (1)a scoping calculation is used to assess the magnitude of the response change; (2)then extra calculations are made which satisfy precision requirements; (3)finally a polynomial fit (linear or parabolic) is constructed from the pool of available calculations and used to estimate the sensitivity coefficient.

Numerous schemes have been designed within this general framework, for example, using one or two scoping calculations and one or two extra calculations, for a range of two to four calculations per input parameter. Clearly with nuclear data one cannot hope to perform DP on all 80,000 parameters. However, CASMO-5MX/DP serves numerous purposes:(1)provide sensitivity profiles for code-to-code comparisons (e.g., with SCALE6 TSUNAMI),(2)provide reference local, first-order uncertainty results to assess other CASMO-5MX methodologies, such as SS,(3)provide sensitivity coefficients for nonnuclear data parameters, for example, fuel enrichment.

Figure 2 shows a flow chart for the CASMO-5MX/DP technique. The basic sequence is to begin with a perturbation factor of unity, that is, 𝑝=1, and perform the nominal calculation. After the base calculation, depending on the specific DP mode chosen (see Table 1), the DP driver will select and perform additional perturbed cases. Using the resultant 𝑦^′(1) from the first perturbed case, the DP driver can now calculate a sensitivity coefficient, S. In the “2-point simple” mode, DP would stop here, using simple finite differences (i.e., linear fit) for the estimate of 𝑆. In the “3-point adaptive” mode, a second calculation is performed with 𝑝 estimated such that the new 𝑦^′(2) satisfies a “small but not too small” criterion; for example, only the three least significant digits show variation. 𝑆 is updated using a linear fit of 𝑦 and 𝑦^′(2). In the “4-point adaptive” mode, one additional perturbed case allows a parabolic fit with 𝑆 estimated as the slope of the fit at 𝑝=1.

Table 1 

Summary of DP modes.

Figure 2 

Flowchart of CASMO-5MX/DP direct perturbation methodology.

Although Figure 2 is shown assuming a single output, CASMO-5MX/DP can effectively produce sensitivity coefficients for all outputs simultaneously, especially with the 4-point adaptive scheme. Figure 2 also makes the distinction that nuclear data perturbations are based on relative perturbation 𝑝 and affect the XS perturbation file, whereas perturbations of general input file parameters result in replacement of 𝑥 in the standard input file with 𝑥^′. Once sensitivity coefficients 𝑆 are available, UQ may be performed using standard first-order uncertainty propagation via (3).

2.3. Stochastic Sampling with CASMO-5MX

The CASMO-5MX stochastic sampling (SS) methodology from [9], shown in Figure 3, uses a very similar framework to the DP methodology (Figure 2). The major differences are summarized below.(1)DP varies a single input parameter at a time (𝑥₀→𝑥^′) whereas SS varies them all simultaneously (𝐱_𝟎→𝐱^′). (2)DP is first a sensitivity analysis technique and with UQ possible through local and first-order uncertainty propagation, whereas SS is first a UQ technique (global and all-order) with approximation due to a finite sample size.(3)Due to the adaptive nature, the robust DP presented requires serial execution of up to 4 cases (although sensitivities of different inputs may be investigated simultaneously) whereas SS is inherently parallel.

Figure 3 

Flowchart of CASMO-5MX/SS stochastic sampling methodology.

The basic sequence in SS (refer to Figure 3) is as follows.(1)Each input is sampled 𝑁 times according to their underlying probability distributions and respecting correlation to other inputs, if any. The 𝑛th sample input set contains is denoted 𝐱^′(𝐧), and note that the main data of the XS perturbation file is just the relative perturbations 𝐩^(𝐧).(2)CASMO-5MX is run 𝑁 times with each set of data; that is, 𝑛=1,2,…𝑁.(3)The distribution of the 𝑁 sets of output is analyzed statistically, for example, with the sample variance.

Note that, in Figure 3, stages of the calculation which result in 𝑁 sets of data/files are shown with a “shadow.”

2.4. Stochastic Sampling with MCNPX/NUSS

In parallel to CASMO-5MX/SS, activities to implement SS in the Monte Carlo code MCNPX have led to the development of MCNPX plus nuclear data uncertainty with stochastic sampling, MCNPX/NUSS, which functions very similarly to CASMO-5MX/SS, except that due to the open nature of the MCNPX ACE library format, it is possible to create perturbed nuclear data libraries and source code modification are not necessary, as shown in Figure 4. As in CASMO-5MX, the same simple random sampling procedure is used but a new tool is needed to apply perturbations 𝐩^(𝐧) to create the perturbed ACE library from the nominal one. Note that the decision to perturb data at the ACE library stage, instead of upstream when data is in the ENDF format, is mainly motivated by the relative ease of access to data in the ACE format. Future versions of MCNXP/NUSS may modify data at the ENDF stage.

Figure 4 

Flowchart of nuclear data UQ with MCNPX/NUSS.

Because the currently used VCM library is based on the SCALE 44-group structure, data perturbations 𝐩^(𝐧) are provided in this structure; however, the system is not restricted to any particular group structure for perturbations. In the “library rewriting” stage, a constant perturbation is applied to pointwise data: 𝑥^′(𝐸)=𝑝_𝑔𝑥₀(𝐸)for𝐸^𝑔≤𝐸<𝐸^𝑔−1, (12) for the perturbation 𝑝_𝑔 of group 𝑔 which ranges from lower to upper energies, 𝐸^𝑔 and 𝐸^𝑔−1. Note that with perturbation of partial cross-sections in the ACE library, the total and absorption cross-section must also be adjusted to preserve consistency in the nuclear data files. The final procedure of the MCNPX/NUSS tool is to systematically supply MCNPX calculations with the generated random ACE-formatted nuclear data files. The MCNPX outputs of interest can be analyzed by the same statistical means as in CASMO-5MX/SS, except for a statistical error term which is inherent to the Monte Carlo calculations. When the distribution of an MCNPX output is characterized, it is important to separate the statistical variance from the variance due to data variations: 𝑣=𝑣_stats+𝑣_data. (13) The magnitude of 𝑣_statsis related to the number of neutron histories in the Monte Carlo calculations and has been estimated to be small compared to data contribution (i.e., nuclear data) for all cases considered here.

3. Results

An overview of the UAM Phase I cases analyzed in this paper is provided in Table 2. Notably, there is no depletion and no soluble boron for any of these cases. In the cell cases of exercise I-1, there is no thermal expansion; however, in the PWR lattice case of exercise I-2, thermal expansion has been assumed which decreases the density and increases the size of all materials. The operating conditions of hot zero power (HZP) and hot full power (HFP) dictate the fuel temperature (𝑇_fuel), moderator temperature (𝜌), void fraction (𝑉), and control rod insertion.

Table 2 

Overview of UAM Phase 1 cases analyzed with CASMO-5MX (C5) and MCNPX/NUSS (MC).

3.1. CASMO-5MX Results

This section presents CASMO-5MX results for both exercises I-1 and I-2. All uncertainty results are in terms of relative standard deviation in percent. For both CASMO-5MX/DP and SS, perturbations are made in the 19-group CASMO-5M group structure, unless otherwise noted. The number of samples used was 𝑁=1000 in all cases. With CASMO-5MX calculations, uncertainty was assumed for all nuclides present in each problem and all reactions available in the SCALE6 VCM library.

3.1.1. Exercise I-1: Cell Physics

The uncertainty summary of exercise I-1 cases is given in Table 3 for the PB-2 (BWR) cases, including results for both CASMO-5MX/DP (C5MX/DP) and CASMO-5MX/SS (C5MX/SS), and in Table 4 for the TMI-1 and Generation III (Gen-III) MOX cases, only with CASMO-5MX/SS. Results show the general trend in eigenvalue uncertainty of approximately 0.5% and 1-group cross section uncertainty of about 1% for most absorption cross sections and nuclides which have mainly thermal fission, but about 4% for nuclides which have significant fast fission. Making the spectrum harder, by introduction of 40% void in the PB-2 HFP case or by using MOX fuel (in the Gen-III case), increases the influence of the fast spectrum, which almost always has higher uncertainty than the nuclear data in the thermal range.

Table 3 

Uncertainty summary for exercise I-1 PB-2 cases.

Table 4 

Uncertainty summary for exercise I-1 TMI-1 and Gen-III MOX cases.

To assess the effect of the perturbation group structure, two additional group structures were investigated as shown in Table 5: the next finer 31-group structure in CASMO-5 M and the 44-group structure of SCALE6. CASMO-5MX/DP was used in order to investigate the breakdown of the uncertainty, that is, which uncertain nuclear data contributed most to an uncertain output. This is presented in terms of the variance fraction, that is, the variance due to that parameter divided by the total variance, which naturally sums to unity.

Table 5 

Effect of changing the perturbation group structure for the PB-2 HZP cell case.

The most influential parameters are easily defined by sorting from greatest to least variance fraction, and the cumulative value can be used to limit the important parameters, for example, the set representing 99% of the total variance, as shown in Figure 5 for the eigenvalue uncertainty and in Figure 6 for the 1-group U-235 fission and U²³⁸ absorption cross section uncertainty. Good agreement between the uncertainty breakdowns is observed except for the U-235 fission spectrum (MT=1018) component which increases considerably with the 44-group structure. As shown in Figure 7, this was found due to the coarse fast groups in the CASMO-5 M 19- and 31-group structures and the highly varying U-235 fission spectrum uncertainty in the fast range in the native 44-groups. Because all perturbations are applied to the CASMO-5MX 586-group library structure, detailed sensitivity profiles may be generated as shown in Figure 7.

Figure 5 

Breakdown of eigenvalue uncertainty as function of perturbation group structure for the PB-2 HZP cell case.

(a)

(b)

Figure 6 

Breakdown of 1-group U-238 absorption cross section uncertainty (a) and U-235 fission cross section uncertainty (b) in terms of variance fractions for the PB-2 HZP cell case.

(a)

(b)

Figure 7 

SCALE6 nuclear data uncertainty for U-235 fission spectrum, MT=1018, (a) and sensitivity profile of the 1-group U-235 fission cross section to the U-235 fission spectrum (b) for the PB-2 HZP cell case.

3.1.2. Exercise I-2: Lattice Physics

The lattice physics cases in exercise I-2 are concerned with propagating both nuclear data uncertainty and the so-called “technological parameter” uncertainty to the two-group nodal data used in conventional core simulators based on two-group nodal diffusion. The output parameters of interest here are mainly the homogenized macroscopic cross sections for fast and thermal absorption (Σ¹_𝑎 and Σ²_𝑎), neutron production (𝜈Σ¹_𝑓 and 𝜈Σ²_𝑓), removal (Σ_rem), diffusion coefficients (𝐷¹ and 𝐷²), and assembly discontinuity factors (ADF¹ and ADF²). A summary of nodal parameters’ nominal values and uncertainties are shown in Table 6 for the TMI-1 PWR assembly at HFP conditions only, with control rods out (unrodded) and in (rodded), considering only nuclear data uncertainty. Additionally, the uncertainty in pin powers was examined at 3 locations: (i) the location of the unrodded case peak power (unr. peak loc.), (ii) the location of the rodded case peak power (rod. peak loc.), and (iii) the gadolinium pin power (Gd pin loc.). See Figure 8 for the locations in the southeast quarter of the 17×17 PWR assembly. The uncertainty in unrodded assembly pin powers was remarkably low; only for the gadolinium pin is the uncertainty greater than 0.1%. In the rodded assembly, pin power uncertainty was slightly greater, on the order of 0.2% for most pins.

Table 6 

Uncertainty summary for the exercise I-2 TMI-1 HFP case, assuming only nuclear data uncertainty.

Figure 8 

PWR assembly locations (southeast quarter shown).

At the time of this publication, the probability distributions of the technological parameters were not generally agreed upon, and so only a sensitivity analysis has been performed using CASMO-5MX/DP which can easily compute sensitivity coefficients of any input file parameter. As in the benchmark specification, five technological parameters were considered: fuel density (fdens), fuel enrichment (enr), fuel pellet radius (rfuel), clad thickness (tclad), and gap thickness (tgap). The sensitivity coefficients with respect to each technological parameter are shown in Table 7. One generally sees the highest sensitivity to the radius of the fuel pellet (rfuel) which can be over 1% variation in an output per 1% variation in pellet radius.

Table 7 

Sensitivity summary with respect to technological parameters for exercise I-2 TMI-1 HFP (values > 0.5 shown in bold).

3.2. MCNPX Results

Results obtained with MCNPX/NUSS for eigenvalue uncertainty in the PB-2 and TMI-1 cell models at HZP are shown in Table 8. Simultaneous variations were performed for U-235 and U-238 reactions, consistent with CASMO-5MX/SS with one exception; the (𝑛,𝛾) partial cross-section (MT=102) is considered explicitly in MCNPX and not the total capture (MT=101) as in CASMO-5 M. For some nuclides with significant (𝑛,𝛼) reactions, comparisons would not be consistent as MT=101 includes (𝑛,𝛼) but MT=102 does not, but for U-235 and U-238, the difference between (𝑛,𝛾) and total capture is minor. Due to the long runtimes of MCNPX calculations, only 𝑁=80 samples were made; however, this achieved statistical uncertainty more than two orders of magnitude less than the data uncertainty for these cases.

Table 8 

Uncertainty summary using MCNPX/NUSS for HZP cases.

Although the number of samples was fairly small at 80, a study of the running average eigenvalue and uncertainty (one-sigma error bars) in Figure 9 shows little sample bias in the sample mean and stable behavior of the sample standard deviation. Additional discussion may be found in [13].

(a)

(b)

Figure 9 

Cumulative moving average of eigenvalue with MCNPX/NUSS versus number of samples for exercise I-1 PB-2 HZP (a) and TMI-1 HZP cases (b).

4. Discussion

In this section, various results from the previous section will be further discussed, namely, (i)BWR uncertainties predicted by both the CASMO-5MX/DP and SS methodologies,(ii)BWR versus PWR uncertainties,(iii)UO2 versus MOX uncertainties,(iv)CASMO-5MX versus MCNPX/NUSS results.

4.1. Comparison of BWR Uncertainties versus UQ Methodology

Consistent trends are observed with both methodologies for the exercise I-1 PB-2 (BWR) case, with slightly higher uncertainties observed at HFP, both in eigenvalue (denoted “K-inf”) and 1-group cross sections, especially U-238 fission. This is due to spectrum hardening in the HFP case, with nearly 40% void, which acts to increase uncertainty because data in the fast range is generally more uncertain. For the 1-group cross sections, a faster spectrum also increases the impact of U-238 inelastic scattering, which contributes greatly to the overall uncertainty [9]. Assuming DP as a reference solution, SS shows excellent agreement (see Figure 10), with smaller eigenvalue uncertainty by less than 4% and larger cross section uncertainty by at most 3% (U-238 absorption at HFP).

(a)

(b)

Figure 10 

Comparison of CASMO-5MX/DP and/SS methods.

4.2. Comparison of Stochastic Sampling Uncertainties versus LWR Reactor Type

At HZP conditions, almost identical uncertainties are observed for the exercise I-1 PB-2 (BWR) and TMI-1 (PWR) cases (see Figure 11). There is slightly larger uncertainty at HZP for the U-235 1-group cross sections due to higher enrichment in the PWR (4.85 wt%) compared to the BWR (2.93 wt%). Because of the previously discussed hardening of the spectrum for the BWR case at HFP, the uncertainty in the U-238 1-group fission cross section is noticeably higher.

(a)

(b)

Figure 11 

Comparison of PWR and BWR uncertainties with CASMO-5MX/SS.

4.3. Comparison of Uncertainties for UO2 and MOX Fuel Types

For MOX fuel from the exercise I-1 Gen-III MOX case, nearly double the uncertainty (0.95%) in eigenvalue is observed compared to UOX fuel (0.51%). See the graphical summary in Figure 12. This marked increase is not only due to the higher uncertainty for the Pu isotopes but also due to the faster spectrum in those cases, which increases uncertainty due to the shift to the more uncertain fast range. Notably, Pu²⁴² and Am²⁴¹ 1-group cross section uncertainties are greater than 4%.

(a)

(b)

Figure 12 

Comparison of MOX and UO2 fuel uncertainty.

4.4. Comparison of Uncertainties for CASMO-5MX/SS and MCNPX/NUSS

The MCNPX results showed a total uncertainty in eigenvalue of 0.54% using MCNPX/NUSS which was very consistent with both the CASMO-5MX/SS and CASMO-5MX/DP results using the same nuclear data uncertainty but different nuclear data libraries and codes systems. Additional tests cases with one-at-a-time perturbations of single reactions have been prepared for a more detailed investigation, comparing to a breakdown from CASMO-5MX/DP, with perturbations in the SCALE6 44-group for maximum consistency with MCNXP/NUSS. The results in Table 9 show the top 5 contributors according to each methodology, and in general one sees excellent agreement. It is perhaps only interesting that CASMO-5MX/DP shows 0.53% uncertainty in the top 5 whereas MCNPX/NUSS shows 0.50%.

Table 9 

Comparison of top 5 contributors for MCNPX/NUSS versus CASMO-5MX/DP for exercise I-1 PB-2 HZP.

5. Conclusions

The UAM benchmark has provided the opportunity to develop state-of-the-art methodologies for uncertainty quantification (UQ) and the framework for international collaboration and comparison. At PSI, within the STARS project, the first development was CASMO-5MX, a modification of the production CASMO-5 M code to perturb nuclear data libraries through an auxiliary input file with the capabilities to provide perturbations in any group structure and perturb individually the inelastic (MT=4) and elastic (MT=2) scattering components despite the internal use of a combined scattering matrix with elastic and inelastic scattering lumped together. Building on this capability, a sensitivity analysis (SA) tool using direct perturbation (DP) was developed, CASMO-5MX/DP, which performs adaptive perturbations in order to robustly estimate sensitivity coefficients of arbitrary outputs with respect to arbitrary inputs, including nuclear data. Using standard first-order uncertainty propagation, CASMO-5MX/DP can also be used for local, first-order UQ. However, to be used for production UQ, CASMO-5MX/DP requires too many calculations, and for these reasons, a second UQ methodology based on stochastic sampling (SS) was developed, CASMO-5MX/SS, which can provide uncertainty estimates for arbitrary outputs at a fixed cost of 100 to 1000 calculations. Most recently, development of an SS methodology for a continuous-energy, Monte Carlo code, MCNPX, was initiated, called MCNPX/NUSS.

Results for the UAM benchmark exercises were presented, including the LWR cell cases from exercise I-1 and the PWR assembly case from exercise I-2. For the cell cases, uncertainty in the eigenvalue and 1-group collapsed microscopic cross sections (in terms of relative standard deviation) was found to be about 0.5% and 1%, respectively. For the Gen-III MOX case, the eigenvalue uncertainty was nearly double (1%) and Pu-242 and Am-241 1-group cross sections uncertainties’ reached 5%. In UO₂ fuel, the most important contributors to the eigenvalue uncertainty were found to be U²³⁸ capture (MT=101), U²³⁵ neutrons per fission (MT=452), and U-235 capture (MT=101), accounting for over 80% of the total variance in eigenvalue. For 1-group cross section uncertainty, U-238 inelastic scattering (MT=4) accounted for well over 50% of variance alone and usually more. For the TMI-1 PWR assembly case, uncertainty in eigenvalue was consistent with the cell cases at about 0.5%. Uncertainty in other assembly outputs ranged from less than 0.1% for the assembly discontinuity factors (ADFs), powers at the nominal peak pin locations, and the thermal diffusion coefficient (𝐷²) to 1% for the fast diffusion coefficient (𝐷¹) and the fast absorption cross section (Σ¹_𝑎). Both rodded and unrodded cases were analyzed and uncertainty was found to remain the same or slightly increase when control rods were inserted.

Finally, sensitivity coefficients were calculated for technological parameters for the exercise I-2 TMI-1 PWR assembly and it was found that the radius of the fuel pellet is the most sensitive parameter, having sensitivity coefficients of absolute value from 1 to 2 for many outputs. For example, the sensitivity coefficient of the removal cross section (Σ_rem) with respect to pellet radius is −2.1 for unrodded case, which means that, for a 1% change in pellet radius, Σ_rem will decrease by 2%! It is clear, however, that a better understanding of the distributions of the technological parameters is necessary, and in particular, how the batch-based nature of manufacturing introduces correlations across the fuel pellets, assemblies. For example, should all the fuel pellets in a single assembly be considered to come from the same batch, different batches, or a fixed number of batches? If a fixed number of batches, how is it determined which pellets are from which batch? Answering these questions requires either more knowledge of how a particular fuel assembly was manufactured or simulation of the actual manufacturing processes! Otherwise, conservative, limiting cases must be created, which is in direct opposition to the overarching goal of best estimate analyses with UQ.

Future work in the area of neutronics UQ at PSI includes enhancement of the MCNPX/NUSS continuous-energy Monte Carlo strategy, implementing the capability to perturb fission product yields and decay constants, and extension of the SS methodology from the lattice code CASMO-5 M to the core simulator SIMULATE-3.