Multiphysics Multiscale Coupling Modeling for Nuclear Reactor and Its Uncertainty Quantification
View this Special IssueResearch Article  Open Access
Jason Hou, Maria Avramova, Kostadin Ivanov, "BestEstimate Plus Uncertainty Framework for Multiscale, Multiphysics Light Water Reactor Core Analysis", Science and Technology of Nuclear Installations, vol. 2020, Article ID 7526864, 18 pages, 2020. https://doi.org/10.1155/2020/7526864
BestEstimate Plus Uncertainty Framework for Multiscale, Multiphysics Light Water Reactor Core Analysis
Abstract
Tremendous work has been done in the Light Water Reactor (LWR) Modelling and Simulation (M&S) uncertainty quantification (UQ) within the framework of the Organization for Economic Cooperation and Development (OECD)/Nuclear Energy Agency (NEA) LWR Uncertainty Analysis in Modelling (UAM) benchmark, which aims to investigate the uncertainty propagation in all M&S stages of the LWRs and to guide uncertainty and sensitivity analysis methodology development. The BestEstimate Plus Uncertainty (BEPU) methodologies have been developed and implemented within the framework of the LWR UAM benchmark to provide a realistic predictive simulation capability without compromising the safety margins. This paper describes the current status of the methodological development, assessment, and integration of the BEPU methodology to facilitate the multiscale, multiphysics LWR core analysis. The comparative analysis of the results in the standalone multiscale neutronics phase (Phase I) is first reported for understanding the general trend of the uncertainty of core parameters due to the nuclear data uncertainty. It was found that the predicted uncertainty of the system eigenvalue is highly dependent on the choice of the covariance libraries used in the UQ process and is less sensitive to the solution method, nuclear data library, and UQ method. HightoLow (Hi2Lo) model information approaches for multiscale M&S are introduced for the core single physics phase (Phase II). In this phase, the other physics (fuel and moderator), providing feedback to neutronics M&S in a LWR core, and timedependent phenomena are considered. Phase II is focused on uncertainty propagation in single physics models which are components of the LWR core coupled multiphysics calculations. The paper discusses the link and interactions between Phase II to the multiphysics core and system phase (Phase III), that is, the link between uncertainty propagation in single physics on local scale and multiphysics uncertainty propagation on the core scale. Particularly, the consistency in uncertainty assessment between higherfidelity models implemented in fuel performance codes and the rather simplified models implemented in thermalhydraulics codes, to be used for coupling with neutronics in Phase III is presented. Similarly, the uncertainty quantification on thermalhydraulic models is established on a relatively small scale, while these results will be used in Phase III at the core scale, sometimes with different codes or models. Lastly, the uptodate UQ method for the coupled multiphysics core calculation in Phase III is presented, focusing on the core equilibrium cycle depletion calculation with associated uncertainties.
1. Introduction
The sensitivity analysis (SA) and uncertainty quantification (UQ) have been integrated into the nuclear reactor design and safety analysis to establish the accuracy and confidence bounds for bestestimate modelling and simulation codes. Traditional approach for nuclear reactor safety involved deliberately pessimistic hypotheses in computations to satisfy required safety margins; however, this approach usually tends to produce excessive conservatisms. In recent years, the accuracy of the computational codes has been improving with the use of more realistic models and hypotheses as well as increased computation resources. The bestestimate approach, typically accompanied with confidence bounds obtained from uncertainty analysis, is regarded in both nuclear industries and regulation as an acceptable alternative to the traditional conservative regulation approach. The Organization for Economic Cooperation and Development (OECD)/Nuclear Energy Agency (NEA) benchmark for Uncertainty Analysis in Modelling (UAM) of Light Water Reactors (LWRs) (known as the LWR UAM benchmark) has been established for over a decade to facilitate the development, assessment, integration, and validation of available SA and UQ methods for bestestimate LWR design and safety calculations, which are currently used in the nuclear power generation industry and regulation. Three main domains of nuclear reactor engineering, corresponding to the three major physics phenomena in reactor systems, are considered in the benchmark, namely, the neutronics, thermalhydraulics, and fuel thermal/mechanical behaviour.
The benchmark first defines several exercises by subdividing the complex LWR systems into smaller problems with varying scales (corresponding to a pin cell, assembly, and core geometry model), each of which will contribute to the total uncertainty of the final coupled system calculation with assumptions made at each geometric level well identified. The uncertainties of the resulting LWR system calculations are determined for each step and are then propagated through an integral LWR system Modelling and Simulation (M&S) for which high quality plant experimental data exist. The benchmark is being carried out in three phases, namely, the neutronics, core, and system phases [1–3]. Three major LWR types are selected based on previous benchmark experience and available data: Pressurized Water Reactor (PWR), Boiling Water Reactor (BWR), and Russian design of PWR (VVER). For each exercise, several test problems are considered in order to cover the different situation targets: reactor types, normal operation conditions, and accidental transients. Significant amount of research and development has been accomplished since the beginning of the LWR UAM benchmark. This paper summarizes the current status of the method development and implementation with a focus on the PWR test cases due to the page limit.
The rest of the paper is organized as the following: Section 2 provides an overview of the two classes of UQ methods that have been developed and applied to the benchmark exercises, as well as the specification of TMI1 reactor related test cases. Section 3 summarizes one of the concluding activities of the standalone neutronics phase (Phase I), that is, the comparative analysis of the results submitted for PWR related cases for understanding the general trend of the uncertainty of core parameters due to the nuclear data uncertainty. The UQ methodologies developed for the core single physics phase (Phase II) are discussed in Section 4, primarily focusing on thermalhydraulics and fuel M&S. Special attention is given to the link and interactions between Phase II and the multiphysics core and system phase (Phase III), that is, the link between uncertainty propagation in single physics on local scale and multiphysics uncertainty propagation on the core scale. Section 5 presents the current status of the uncertainty propagation and quantification method development for the multiphysics coupled PWR core simulation. Last, the general conclusion of the study and future work are provided in Section 6.
2. Uncertainty Quantification Methodologies
2.1. Sources of Uncertainties
In principle, the sources of input uncertainties in computer code simulations include data uncertainties (e.g., nuclear data, geometry, and materials (manufacturing or technological uncertainties)); uncertainties in physical models and numerical methods (e.g., approximations in the numerical solution, nodalization, homogenisation approaches); and uncertainties due to the imperfect knowledge of boundary and initial conditions. Input uncertainties are defined for each exercise as a combination of new uncertainties and the ones propagated from previous exercises.
2.2. Uncertainty Quantification Methodologies
In general, two sets of uncertainty quantification (UQ) methods were pursued by the participants in the LWR UAM benchmark, namely, the deterministic and stochastic methods.
The deterministic method calculates the sensitivity of the system response with respect to uncertain input parameter using Perturbation Theory (PT) and computes an estimate for the covariance matrix by linearizing the response . Here, is the response vector sensitivity matrix. With the linearization, the covariance matrix can be calculated by folding sensitivities with the Variance and Covariance Matrix (VCM) of input parameters:where denotes the VCM of the input parameter . Equation (1) is known as the firstorder uncertainty propagation formula or the socalled “Sandwich rule” [4]. The sensitivity matrix can be calculated using various codes and an example is given in [5].
The sampling method relies on the sampling of the uncertain input parameters provided in the VCM and statistically analyses the calculated output responses. The variance is computed usingwhere is the number of samples and is the sample mean of the response. Examples of the implementation of the statistical sampling can be found in [6, 7].
The potential difference between the calculated uncertainty using deterministic and sampling methods can often be justified by the relative confidence intervals, which is related to the sampling approach. For a given sample size, the confidence bounds of uncertainties of core output responses could be quantified if normality assumption is valid. If a simple random sample size N is obtained from a normally distributed population with true mean µ and true standard deviation , then the following constructed parameter follows a chisquare distribution with N − 1 degree of freedom:where is the sample variance. For a twosided uncertain parameter, the criterion of 95% confidence level indicates that 95% of the value is bounded inside the interval of [,]. The confidence interval of the standard deviation can be derived as
Figure 1 shows the upper and lower limit of the 95% and 99% confidence level are depicted as a function of the sample size. For a twosided distribution, the 95% confidence interval of the is evaluated as [, ], which covers the true variance of R with 95% confidence level. Therefore, the 95%/95% confidence interval of the true uncertainty can be mathematically determined as [90% σ_{p.}, 113% σ_{p.}].
2.3. TMI1 Specification
This paper summarizes the recent advancement in the LWR UAM benchmark, focusing on the PWR exercises based on the Three Mile Island Unit 1 (TMI1) reactor, which is used as an example to illustrate the developments and implementations performed with the LWR UAM benchmark framework.
The TMI1 core consists of 177 fuel assemblies, each of which contains 208 fuel rods, 16 guide tubes, and 1 instrumentational tube. There are 11 types of fuel assemblies in the TMI1 active core with various fuel enrichment (4.00%, 4.40%, 4.85%, 4.95%, and 5.00%) and configurations with regard to the configuration of the burnable poison (BP), gadolinium pins (GdO_{2} + UO_{2}), and control rod banks. Figure 2 depicts one octant of the core radial layout, where assembly H8 is located in the core center. Detailed geometry setup, material composition, and enrichment can be obtained from LWR UAM benchmark specifications [1, 3].
3. Uncertainty Propagation in Reactor Physics Simulations
The first phase of the LWR UAM benchmark [1] is dedicated to the standalone multiscale neutronics M&S and carried out in three steps, each corresponding to one of the steps of the standard LWR simulation approach: cell physics (to produce multigroup microscopic cross section libraries), lattice physics (to derive multigroup homogenized macroscopic cross section sets), and core physics (to assess fullcore performance parameters). The sources of the input uncertainties considered in Phase I exercises include the neutron cross section data, supplemented by the VCMs, and asbuilt manufacturing uncertainties in material composition and geometric dimensions.
Most participants of the Phase I exercises focused on the quantification and propagation of nuclear data induced uncertainty; thus the VCM is simply the nuclear data covariance information, available either by processing covariance data files provided by major nuclear data libraries (NDLs) or in the SCALE code package [8]. Two SCALE VCMs have been used: the 44group library distributed with SCALE 6.0 and SCALE 6.1 and the updated library available in SCALE 6.2 which is based on 56group and 252group structures [8]. The former contains uncertainty data for 401 materials with important isotopes taken from highfidelity nuclear data evaluations including ENDF/BVII.0, ENDF/BVI, and JENDL3.3. The latter is based on the EDNF/BVII.1 data for 187 nuclides, combined with information of ∼215 nuclides from the SCALE 6.1 VCM. Note that in many cases, the SCALE VCMs were transformed into a user specified energy group structure to satisfy the requirement of the neutronics code used in the UQ and propagation process [9]. Table 1 provides a summary of numerical calculations carried out to obtain the results corresponding to the neutronics cases, including the NDL, transport code (solution method), VCM, and UQ method.
 
^{1}Covariance data were collected from various sources, including JEFF3.2, ENDF/BVII.1, JENDL4.0, and TENDL2011. 
3.1. Cell Physics
For each test case in Exercise I1, the following results are requested from the participants: calculated kinf and its associated uncertainty, the top five neutronnuclide reactions that contribute the most uncertainty to kinf, and covariances of selected onegroup cross sections generated in the pin cell calculation.
A twodimensional (2D) model adopted from the TMI1 reactor was chosen as the representative PWR test problem. It is fuelled with 4.85% enriched UO_{2} fuel, and both Hot Zero Power (HZP) and Hot Full Power (HFP) conditions were specified. In total, 22 sets of results have been submitted for the HZP conditions, from which calculation related information was collected to facilitate the comparative analysis of the results in order to reveal the relation between the choice of calculation parameters (e.g., transport solution method, covariance library, and UQ method) and output uncertainties. Figure 3 shows the predicted kinf and associated uncertainties. The nominal value of kinf spans a range of ∼1200 pcm, which is believed to be a result of the differences among neutron transport solvers and base nuclear data libraries involved in the calculation and/or potential modelling errors. Further investigation is required to quantify the contributions to this rather large deviation. The averaged nominal value and Relative Standard Deviation (RSD), or , of all predicted kinf are 1.43019 and 0.52%, respectively.
More detailed analyses were performed to determine the correlation between each of the calculation parameters. The uncertainty of kinf and the results indicate that the choice of the covariance library strongly impacts the RSD of kinf, as shown in Figure 4, while other parameters have limited influence. The average RSD of kinf calculated using the SCALE 6.0/SCALE 6.1 covariance libraries is below 0.5% (0.47% and 0.48%, resp.), while the value corresponding to SCALE 6.2 is 0.54%. Only one dataset each was submitted using the SCALE 5.1 and JENDL4.0 covariance libraries and the RSDs are ∼0.5%. The ENDF/BVII.1 covariance data yields the highest kinf RSD of 0.61%, despite the fact that the one dataset with extreme low uncertainty is included. It is difficult to determine the reason of low uncertainties presented in the comparison (e.g., Cases 7 and 13) based on the submitted data, and hence further investigation will be performed.
The trend above can be explained by the difference of the covariance information in various sources. On one hand, the covariance data in the ENDF evaluations are generated as part of the cross section evaluation process and represent uncertainties and correlations in differential data. The use of this covariance to calculate uncertainties for integral quantities such as kinf will usually result in an overestimation of the uncertainty. On the other hand, a range of different tests, such as the critical benchmark experiments, are performed to investigate and adjust the new covariance data from the NDL before being incorporated into the SCALE covariance libraries [9]. For example, two of the changes that have significant impacts on many experiments are the modifications made to the covariance data for nubar (average number of neutrons per fission reaction ) of ^{235}U and ^{239}Pu, with the first one relevant to the case study shown in this paper. The uncertainty of ^{235}U nubar in the thermal range increases from 0.31% in SCALE 6.1 to 0.7% in ENDF/BVII.1, which is responsible for differences exhibited in the covariance testing in the lowenriched, water moderated uranium oxide pin array (LCT) systems. Consequently, it was reduced to 0.39% in SCALE 6.2 library, which is consistent with the value in JENDL3.3 [8]. This explains why the calculated RSD of kinf using SCALE 6.2 data is slightly higher than that using SCALE 6.0/SCALE 6.1, while the value corresponding to ENDF/BVII.1 is the largest.
It is convenient to use the PT method to compute the sensitivity coefficients of output variables with respect to nuclear data as compared with the sampling approach, thus making it possible to determine the most influential nuclidereaction pair to the predicted kinf uncertainties by sorting them from greatest to least variance fraction. For the TMI1 HZP case, 13 sets of submitted results include such information and Figure 5 shows the occurrence of various nuclidereaction pairs as the top 5 contributors. Although a certain degree of diversity can be found in the ranking as it includes up to 10 nuclidereaction pairs, some reactions dominate the contribution to the uncertainty of kinf, such as ^{238}U capture, ^{235}U nubar, and ^{235}U capture.
By definition, these main contributors to the uncertainty are identifiable due to (1) the highest sensitivities associated with such reactions, or (2) the highest value of the associated covariances, or (3) a combination of both. For example, the kinf is quite sensitive to the ^{238}U capture cross section, especially in the unresolved resonance regions, while the evaluated cross sections exhibit large uncertainties [10]. This is why, on one hand, ^{238}U capture reaction is the predominant component to the total uncertainty when covariance libraries of SCALE 5.1/SCALE 6.0/SCALE 6.1 are utilized. On the other hand, ^{235}U nubar tops the ranking if the SCALE 6.2 or ENDF/BVII.1 covariance libraries are used, which again confirms the evolution of the ^{235}U nubar uncertainty as described above. The ^{235}U capture is another important contributor as it almost always ranks the third regardless which covariance library is used. The sensitivity profile of the three most influential nuclidereaction pairs is depicted in Figure 6. It is worth mentioning that some nuclidereaction pairs only associate with the ENDF/BVII.1 covariance library, including the ^{1}H capture, ^{238}U nubar, and ^{238}U elastic scattering, which can also be explained by the difference in covariance libraries. For example, the ^{1}H thermal capture in SCALE covariance library is adopted from JENDL 3.3 [11], which is lower than that in ENDF/BVII.1 by a factor of five.
3.2. Lattice Physics
The main focus of this exercise is placed on the propagation of nuclear data uncertainties through lattice calculations to the uncertainty of target output variables, primarily the fewgroup constants (homogenized cross sections and other nodal parameters). The other sources of uncertainty to be considered in this exercise are the uncertainties associated with methods and modelling approximations embedded in the lattice codes.
The PWR lattice model is a 2D fuel assembly of 15 × 15 rods with 1 central instrumental tube, 16 guide tubes, and 4 corner Godlinia pins containing integral Gd burnable poisons. The case was modelled under both unrodded and rodded conditions at HZP and HFP. There are 12 and 11 sets of results submitted for the HZP unrodded and rodded case, respectively.
Figure 7 compares the predicted kinf uncertainty of the unrodded and rodded lattice, from which it can be seen that the RSD in the unrodded case is slightly lower than that in the rodded case, regardless of the choice of the VCM, mainly due to the higher nominal value of kinf when the absorber material is absent in the unrodded model. A similar trend of the VCM’s influence is observed as in cell physics results; that is, the use of SCALE 6.2 leads to larger predicted uncertainties (after removing one obvious outlier from the results), which is primarily caused by the large uncertainty of the ^{235}U nubar.
(a)
(b)
The uncertainties of the homogenized twogroup nufission cross sections are compared based on the VCM, as shown in Figure 8. Again, the higher nuclear data uncertainty of the nubar reaction in SCALE 6.2 is responsible for the large value of RSD calculated by Case 23. It can also be found that the uncertainty of homogenized nufission cross section is more profound in the fast group than the thermal group consistently. For example, the value is ∼0.62% in the fast group versus ∼0.49% in the thermal group in Case 23. It is known that the ^{235}U nubar covariance is lower in the fast and epithermal range than in the thermal energy range in the SCALE 6.2 VCM [8], so are the sensitivity coefficients of system kinf with respect to the ^{235}U nubar and fission cross section [12]. The comparison reflects the fact that the nominal value of nufission in the thermal group (0.19 cm^{−1}) is more than 10 times larger than that in the fast (9.1 × 10^{−3} cm^{−1}) group.
(a)
(b)
The relationship between the kinf, flux, and group constants (cross sections, diffusion coefficients, and assembly discontinuity factors) in the twogroup representation can be determined by calculating the correlation coefficients between variables. Selected correlation coefficient matrices are depicted in Figure 9, in which red and blue colours represent positive and negative correlation, respectively, between two parameters, while correlation close to zero, in white, indicates that there is a low correlation.
(a)
(b)
(c)
(d)
On one hand, the matrices agree well with each other and similar trend is observed. For example, the fast flux is negatively correlated, and the thermal flux is positively correlated, with Group 1 absorption, scattering, and transport cross sections. A higher positive correlation is found between kinf and Group 2 nufission cross section. On the other hand, some unique features are also observed only in one of the cases. Case 12 shows strong crosscorrelations between the scattering and transport cross sections of both thermal and fast group. As a result, large correlation coefficients are shown for the total cross section in Groups 1 and 2. In Case 6, the values of kinf related correlation coefficients are generally larger for cross sections and fluxes, as compared with the other three cases. Further investigations are necessary to identify the reason of those features, that is, whether it is associated with the VCM used in the calculation; however, it is expected that the difference will be propagated with the group constants to the core calculation and will affect the uncertainty of the core parameters.
3.3. Core Physics
The uncertainties of the fewgroup cross sections evaluated in Exercise I2 are to be propagated through the standalone neutronics core calculation in Exercise I3 to the parameters of interest such as the core keff and power distribution. Various uncertainty propagation methodologies can be used, as summarized in Table 2 [12]. In the full deterministic approach, the PT calculation is performed at both lattice and core level, and the VCM of the fewgroup homogenized constants generated in the lattice calculation is used to evaluate the uncertainty of core responses. The onestep approach relies on the stochastic sampling method on both lattice and core levels and is named as such because it involves a onetoone connection between lattice calculations to generate the fewgroup cross section library and the core calculation that read this library as input [13]. The twostep approach combines the stochastic sampling and deterministic methods. Detailed explanation and a complete list of references can be found in [12].

It should be noted that all the aforementioned approaches follow the standard reactor simulation procedure, which involves the generation of homogenized constants, simplification of core geometry, and application of lower order solvers such as the nodal diffusion method. The major modelling difference occurs on whether the spatial homogenisation is performed over the assembly or pin cell. Either way, this procedure will inevitably introduce additional discrepancies to the calculated uncertainty of core responses.
The PWR model defined in the core physics is the TMI1 fresh core at HZP state with all control rods fully inserted. In total, 9 sets of results have been collected. Most results were obtained using the onestep and twostep approaches mentioned above, such as in [14], except for Cases 11 and 35, where the direct fullcore transport calculation was performed with explicit modelling of the geometry. The estimated relative uncertainties of the core eigenvalue are shown in Figure 10, and the average value over all 9 sets of results is 0.49%, which is close to the observed values in the pin cell and lattice calculations. Also similar to the previous cases is the dependency of the uncertainty on the VCM choice; that is, the calculation using the SCALE 6.2 VCM tends to yield higher uncertainties compared with other libraries, primarily due to the large uncertainty in the neutron production rate. It should be noted that the input uncertainty in the fullcore calculation has a wider spectrum than in the previous steps because of more complicated material compositions (e.g., addition of neutron absorber) and increased approximation to the neutron transport solution.
4. Uncertainty Propagation in Fuel Behaviour and ThermalHydraulics Simulations
4.1. Interactions between Multiscale Single Physics and Multiphysics Modelling
In Phase II, the other physics (fuel and moderator), providing feedback in a LWR core, and timedependent phenomena are considered. Phase II is focused on uncertainty propagation in single physics models which are components of the LWR core coupled multiphysics calculations. Phase III is focused on propagation of multiple uncertainties in coupled multiphysics steady state, cycle depletion, and transient calculations. The interactions between Phase II and Phase III are shown in Figure 11.
Exercise II2 considers timedependent neutronics phenomena–depletion (slow time phenomena) and kinetics (fast time phenomena). This exercise develops also methodologies for cross section modelling, that is, for generation of parameterized cross section, kinetics, and other nodal parameters’ uncertainty libraries to be used in multiphysics calculations in Phase III. Such methodology developed using the Sampler/Polaris sequence of SCALE 6.2.1 [8] is shown in Figure 12.
4.2. Fuel Modelling
Exercise II1 is dedicated to fuel physics in steady state and transient conditions. The focus is on consistency in uncertainty assessment between fine models implemented in fuel performance codes and the rather simplified lumped models implemented in thermalhydraulics codes, to be used for coupling with neutronics tools in Phase III [15].
Highfidelity fuel performance codes are used to model a single fuel rod, using detailed technological and operation data (geometry, enrichment, burnup, etc.) and highfidelity models for fission gas release, cladding corrosion, swelling, and so forth. These codes are used in Exercise II1, while lowfidelity models, embedded in thermalhydraulics codes, are usually preferred for fullcore computation in Phase III. In this way Exercise II1 propagates uncertainties in fuel physics consistently using highfidelity fuel performance codes. The focus is on manufacturing, boundary conditions, and subset of modelling (material properties) uncertainties. Figure 13 shows comparisons of benchmark participants’ results using highfidelity fuel performance codes for the PWR depletion test case based on TMI1 fuel rod design. Exercise II1 includes also a special test case (modelling of one axial node/rodlet of single pin) to evaluate the capability of simplified lumped fuel rod models of system and subchannel thermalhydraulics codes to predict fuel temperature as compared to fuel performance codes.
The socalled Doppler temperature is the main fuel physical feedback quantity used in fullcore coupled multiphysics calculations. The Doppler temperature is a kind of “fuel equivalent temperature” which is used to calculate the Doppler feedback. Several methods exist to calculate it from the fuel pellet temperature distribution. A common method is to apply a weighted average between the centreline temperature and the surface temperature.
Actually the fuel temperature is calculated either for an average lumped fuel rod at assembly (or nodal) level (in the simplified lumped models of thermalhydraulics codes) or for an individual rod in a pinbypin (or cellbycell) calculation (in fuel performance codes). In both cases, the fuel temperature is a result of a onedimension heat conduction problem solution. In Phase III in the coupled multiphysics calculations, the temperature is calculated for each node of the core (3D distribution radial and axial directions) using tabulated fuel properties, the nuclear/thermal power given by neutronics, and the local thermalhydraulic conditions as solved by the core thermalhydraulics model. The rodtocoolant heat transfer coefficient is also given by thermalhydraulics or defined by the user as input parameter. The heat conduction problem is actually solved in two different regions: the cladding (without heat source) and the fuel pellet. The gap is modelled by a thermal resistance (or conductance), directly provided as an input data by the user, or calculated by a dedicated correlation. The fuel and cladding specific thermal properties (density, heat capacity, and conductivity) are also provided, usually as functions of the temperature and of the local burnup.
Since in a highfidelity fuel performance code a single rod is modelled, the fuel temperature is calculated at each axial and radial node of the meshing of the fuel rod, and the resolution of the thermal problem takes into account all the possible physical phenomena that affect the fuel geometry and thermal properties. In summary, in a fuel performance code, fuel temperature field, fuel rod geometry, and gap conductance can be assessed. Uncertainty propagation can be performed on the three most important modelling quantities: fuel thermal conductivity, cladding thermal conductivity, and gap conductance. To summarize, the calculations performed within the framework of Exercise II1 with higherfidelity fuel performance codes help to develop parameterized values and associated uncertainty bounds for these three parameters to be used in the coupled multiphysics calculations in Phase III. The following data for the gap conductance, fuel, and claddings thermal conductivities is generated using highfidelity fuel performance codes—(i) the mean value and (ii) associated uncertainty—both as functions of the following parameters: for cladding conductivity k_{c}(T), as function of temperature (T); for fuel conductivity k_{f}(Bu, T), as function of burnup (Bu) and temperature (T); and for gap conductance, h_{g}(Bu, P) as function of burnup (Bu) and power (P). These data are computed and parameterized for each type of fuel rod, that is, for each type of LWR fuel assembly of interest. We will refer the generation of this data as a preparation of Hi2Lo model information plus uncertainties using highfidelity fuel models for lumped simplified lowerfidelity fuel models. The lookup table representation concept [2] is appropriate for steady state and cycle depletion calculations while for transient/accident simulations the hightolow model calibration of selected coefficients in dynamic gap conductance correlations is preferred [16]. Please keep in mind that the uncertainty values could be extracted also from the experimental datasets and associated uncertainties used to develop the corresponding correlations utilized in the fuel models. As an example, the uncertainty associated with the IAEA correlation for cladding conductivity is shown in Figure 14.
4.3. ThermalHydraulics
As for fuel thermal and mechanical behaviour, there might be different thermalhydraulic models used between Phases II and III. The uncertainty quantification on thermalhydraulic models is established on a relatively small scale, that is, rod bundles, in Exercise II3, while these results will be used in Phase III at the core scale. Detailed twophase flow models, either at subchannel level or at Computational Fluid Dynamics (CFD) level, are used to model a single fuel bundle in Exercise II3, while rather simplified models at fuel assembly level are utilized for fullcore thermalhydraulics in Phase III [15].
Considering that the feedback effect of thermalhydraulics on core behaviour leads to identifying three main parameters of interest: the coolant temperature and/or void fraction and the cladding temperature (i.e., boundary condition for heat conduction problem in fuel rods) resulting from the heat transfer coefficient between the fuel rod and the coolant. These parameters highly depend on the flow: singlephase turbulent flow, nucleate boiling, or film boiling (postCHF conditions).
Hence, Phase II should provide uncertainties on these parameters but also on the transition criteria between the three flow regimes. Due to the large range of available models for twophase flow, the benchmark specification cannot address input uncertainties to be considered for each model and all closure laws. Nevertheless the consistency between Phases II and III should be verified by the participants when different codes or models are used (e.g., fourequation model versus twofluid model). The methodology is the following: the modelling uncertainties have to be tuned in order to have the same effect on the coolant temperature, for example, at subchannel level (for which it is possible to compare the computations to the measurements).
As example, the obtained results of statistical uncertainty analysis performed with subchannel code CTF coupled with the statistical analysis tool, Design Analysis Kit for Optimization and Terascale Applications (DAKOTA) [17] for TMI1 PWR assembly at normal operation conditions as part of Exercise II3, are shown and discussed [18]. The statistical method repeatedly runs the code with inputs sampled from the specified input distributions of uncertainties in boundary conditions, geometry, and modelling. Through statistical analysis of the results, conclusions are drawn about the behaviour of the quantities of interest. For the uncertainty quantification, a Latin Hypercube Sampling (LHS) method with a sample size of 2000 was used. This number of samples was selected by utilizing the maximum proportion of reasonably available computational resources. The TMI1 quarter assembly model at normal operation conditions is utilized. This case is unique because there is not significant void generation in the bundle; therefore, the fluid temperature is selected as a quantity of interest in place of the void fraction. The axial fluid temperature distribution in selected subchannel, axial cladding surface temperature distribution of selected rod, and exit fluid temperature distribution are all shown in Figure 15. An expected result of the singlephase nature of this case is that the standard deviations for all temperatures are approximately constant, especially within each individual axial distribution. The physical modelling of singlephase flow is simpler and better understood than twophase modelling, which leads to smaller overall uncertainty in thermalhydraulic simulations.
(a)
(b)
(c)
In has been shown previously that CFD codes could be used to generate Hi2Lo information with uncertainties for spacer grid and incore structure modelling [19]. These CFDinformed models are developed using a physicsbased approach which utilizes highfidelity data to inform and calibrate a lowerfidelity simulation. It is important to consider only once input uncertainties. For example, the spacer grid effects (position and mixing) are investigated in Exercise II3 for a given type of fuel bundle—in this case, TMI1 PWR bundle. The resulting uncertainty will be considered later in Phase III in the overall uncertainty on investigated output parameters. The methodology utilizes CFDbased data for: griddirected crossflow model; gridenhanced turbulent mixing model; gridenhanced heat transfer model; and spacer grid pressure losses model. This data while utilized in lowerfidelity core thermalhydraulics modelling as part of core multiphysics calculations in Phase III is supplemented with uncertainties propagated through the highfidelity local CFD simulations. These uncertainties are part of input uncertainties propagated through global core multiphysics calculations.
5. Uncertainty Propagation in Multiphysics Reactor Core Simulations
5.1. Approach to Multiphysics Uncertainty Propagation
Phase III of the LWR UAM benchmark is focused on the prediction of key reactor parameters associated with LWR multiphysics simulations. Coupled fuel rod, thermalhydraulics, and neutronics simulations are included by taking into account coupling/feedback effects between the three phenomena. Most of the progresses made so far are on Exercise III1, which aims to identify and propagate input uncertainties through core multiphysics calculations.
Selected output uncertainties from various physics domains/phenomena have been evaluated in Phases I and II of exercises and will be used as input uncertainties in Exercise III1. Those uncertainties can be given in the form of distributions, a single value, or multiple data files and will be combined with other sources of uncertainties in core multiphysics calculations. Following the multiscale, multiphysics modelling strategy defined in the current benchmark, one possible approach for propagating those uncertainties is proposed and demonstrated in Figure 16. The input uncertainty, as discussed below, stems from three physics domains of LWR modelling: neutronics, fuel modelling, and thermalhydraulics.
The neutronics uncertainties are embedded in the twogroup cross section libraries, which include parameterized cross sections, discontinuity factors (DFs), burnup, and kinetics parameters. The fuel modelling uncertainties are applied to the standard lumped fuel rod models of system and subchannel codes via three parameters related to the heat transfer from fuel to coolant: fuel thermal conductivity k_{f}, gap conductance h_{g}, and cladding thermal conductivity k_{c}. In general, the input uncertainties for the core thermalhydraulic model include boundary conditions uncertainties, geometry uncertainties, and modelling uncertainties. At the first stage of implementing the method depicted in Figure 16, the uncertainties associated with different physics spaces are considered independent from each other and hence propagated separately towards the core simulation, although the simulation itself is carried out using the coupled multiphysics code system.
5.2. Uncertainty Propagation for Cycle Depletion Calculations
The TMI1 equilibrium cycle depletion case was specified in the LWR UAM benchmark, Exercise III1, where the reactor core is depleted at HFP of 2771.9 MWth, with the average fuel temperature of 921 K. The mass flow rate through reactor core is assumed to be 16546.04 kg/s with the coolant inlet temperature of 562.67 K. Control rod Groups 1–6 are completed withdrawn, while Group 7 and the axial power shape rod (APSR) Group 8 are 90% and 30% withdrawn, respectively. The core is depleted for 664 EFPD before the critical boron concentration 5 ppm is reached.
The TMI1 core was modelled using PARCS [20] with the fixed mesh size. The core is radially partitioned into 177 nodes based on the onenodeperassembly configuration plus the radial reflector. Axially, the core is discretized into 24 equal nodes in the active core region and 1 node each for top and bottom reflector. The cross sections were generated as a function of fuel temperature, moderator density, boron concentration, control rod insertion, and depletion for 11 types of fuel assemblies. Cross sections for the 3 types of reflectors were assumed to be invariant as those generated under nominal state condition. The twogroup homogenized cross sections were generated based on the ENDF/BVII.1 library using Polaris, which is a 2D LWR lattice physics module in the SCALE 6.2.1 code system [8]. The fluid dynamics and heat transfer calculations were performed using the PATHS (PARCS Advanced ThermalHydraulics Solver) code [21] to provide thermalhydraulics feedback for steady state and depletion calculation via the tight coupling with PARCS. The axial nodalization and spatial mapping (overlays) of thermalhydraulics (TH), neutronics kinetics (NK), and heat structure (HTSTR) models in the coupled depletion simulation is shown in Figure 17.
The nominal depletion calculation was first performed using PARCS/PATHS and the results were compared against the reference solution provided in LWR UAM benchmark specification. With a cycle length of 664 EFPD, critical boron concentration at EOC is 38 ppm, which compares favourably with the reference solution with an error of 33 ppm. Moreover, the maximum relative error of the assembly averaged burnup is −1.39%, as shown in Figure 18.
The uncertainty propagation and quantification were carried out following the approach described in Figure 16. On the neutronics side, this study utilized Sampler (represented by “sampling module 1”) [8], a module for statistical uncertainty analysis in SCALE code package, to sample probability density functions (PDFs) defined in the SCALE 56group covariance library and generate random samples of nuclear data for the subsequent lattice calculation. The fewgroup constants were then generated and converted by the lattice physics code Polaris [8]. On the thermalhydraulics and fuel modelling side, DAKOTA (shown as “sampling module 2”) [22] was employed to produce samples of selected parameters according to their PDFs (see Table 3), using the Latin Hypercube Sampling (LHS) method. The LHS is a stratified sampling technique for which the range of each uncertain variable is divided into segments of equal probability; as a result, it requires less numerous samples than the random sampling method to achieve the same accuracy in statistics. The perturbed input variables for the multiphysics simulator PARCS/PATHS were prepared by pairing one set of fewgroup constants with one sample of fuel modelling and thermalhydraulics parameters. For specific tolerance limits, the number of code runs (or the sample size) can be determined using Wilks’ method [23]:where is the uncertainty, is the statistical confidence, and n represents the sample size. A total number of 146 sets of input variables were formulated in order to ensure that 95% of the results fall in a confidence level of 95%.

A pythonbased interface was developed to streamline the calculation and data processing, including feeding the input parameters, invoking core depletion calculation for 664 EFPD using PARCS/PATHS, and performing the statistical analysis of the collected output responses. The results are shown below.
Figure 19 shows the boron letdown curves of the unperturbed and perturbed depletion calculation as well as the associated statistical results. It can be seen that the peak boron concentration occurs at BOC for all sample runs with the values of 1786 ppm and 1809 ppm for the nominal case and sample mean, respectively. It also shows that the peak critical boron concentration will be bounded within the range of [1656, 1932 ppm] for the tolerance limit of 95% at a confidence level of 95%. Note that the critical boron concentration reduced to 5 ppm before reaching 664 EFPD in some perturbed cases, and the cycle lengths in those cases are thus shorter than others.
The radial and nodal power peaking factors are depicted in Figures 20 and 21, respectively, during the entire cycle depletion. In general, it can be found that the uncertainty of the peak relative power is larger at BOC than EOC. Among all perturbed cases, the 95^{th} percentile of the distribution was calculated at each burnup step for the peak relative power to establish a tolerance limit with 95% confidence. The results are 1.395 and 1.777 for the peak radial and peak nodal power, respectively.
The assembly averaged burnup was also calculated with associated uncertainties, as shown in Figure 22. Since the assembly burnup is not considered as a safety parameter, the 95%/95% confidence interval was evaluated to quantify its uncertainty, as shown in (3), in addition to the sample mean and standard deviation. The maximum uncertainty of the assembly burnup is found to be 0.82% in assembly L13, which is equivalent to 0.84 GWD/MTU, and the 95%/95% confidence interval of the true uncertainty can be estimated as [0.73%, 0.93%].
The mean and upper bound (computed as mean value plus the deviation between the mean and the 95% intervals) of the cycle length were also calculated for all samples, which was determined when critical born concentration reaches 5 ppm, since it is shorter than 664 EFPD in some perturbed cases. By separately perturbing nuclear data and each of the fuel modelling and thermalhydraulics parameters, the total uncertainties observed can be decomposed and the impacts of various uncertainty sources can be evaluated, as shown in Table 4. It is found that the total uncertainty of the cycle length is 36.2 days and mostly contributed by nuclear data uncertainty. Fuel modelling parameters are responsible for uncertainty of only 1.2 days, which is mostly dictated by the uncertainty of the gap conductance. Thermalhydraulic parameters are responsible for uncertainty of 12.3 days, mostly stemming from the core power uncertainty.

The results of the uncertainty breakdown indicate that the uncertainties of the core responses in the cycle depletion calculations are influenced by the input uncertainty of all three physics domains. This is consistent with the conclusion of previous studies, such as [14], that the uncertainty of power and reactivity is mostly dominated by the nuclear data uncertainty, while that of the temperature is strongly dictated by the uncertainty of fuel modelling and thermalhydraulics variables. It is obvious that the evaluation of the input uncertainties is of paramount importance for determining the uncertainty of the responses.
6. Conclusions and Future Work
This paper summarizes the current status and outcome of the development of the BestEstimate Plus Uncertainty (BEPU) framework for the multiscale, multiphysics LWR core analysis under the guidelines of the LWR UAM benchmark, primarily focusing on the PWR cases based on the TMI1 reactor. The standalone neutronics exercises in Phase I are mainly concerned with the propagation of input uncertainties through the standard multistep LWR simulation procedure to key core parameters, such as the multiplication factor. It was observed that the uncertainty estimates for the system eigenvalue due to nuclear data in all scales (the pin cell, lattice, and core) are similar with the Relative Standard Deviation of approximately 0.5% Δk/k. Comparative analyses of all submitted results were then carried out to investigate the dependence of the predicted uncertainties of crucial core parameters on the choice of solution methods, nuclear data libraries, VCMs, etc. It was found that the kinf (keff) uncertainty strongly depends on the VCM and the results generated using the VCM of ENDF/BVII.1 and SCALE 6.2 lead to significantly larger values than those using other VCMs. The sensitivity analysis suggests that the large deviation of the ^{235}U nubar uncertainty in the thermal energy range is responsible for the disagreement. As a followup to the current work, similar analyses will be applied to benchmark cases belonging to other reactor types, including BWR, VVER, and GenIII reactors. It was expected to help identify important nuclide/reaction types to the calculated uncertainty and make recommendations to the further improvement of the nuclear data evaluation for the nuclear reactor application. Other sources of input uncertainties, such as the manufacturing tolerances, will also be taken into consideration in the future comparative analysis.
Using the experience gained on Phase I, a consistent approach is defined in Exercise II2 to parameterize the neutronics uncertainties in a form of fewgroup libraries to be propagated from Phase II to Phase III. Exercises II1 and II3 are used to characterize the uncertainties of highfidelity fuel and thermalhydraulics models regarding parameters of interest at core level, such as nodal Doppler temperature and coolant temperature. These uncertainties are defined as functions of local operating conditions at core level (pressure, burnup, etc.) in order to be propagated in Phase III with rather simplified models for fuel behaviour and core thermalhydraulics. For the three selected fuel physics quantities, hightolowfidelity model information approach is applied, and the obtained uncertainties are propagated in the core multiphysics calculations. Phenomena Identification and Ranking Table (PIRT) analysis will be used for thermalhydraulic quantities where common crosscutting input has the strongest uncertainties impact for the envisioned initial steady state, cycle depletion, and transient applications. It is only in those quantities that the uncertainties will be propagated for fullcore calculation. For such propagation, the experience accumulated in Exercise II3 test problems is utilized. Finally, the multiphysics exercises of Phase III will result in safety and design quantities with the bestestimate value and the overall uncertainty.
In Phase III, the uncertainty quantification of the core equilibrium cycle depletion calculations was investigated. The core depletion calculation was performed using the PARCS/PATHS code, while Sampler and DAKOTA were employed to propagate input uncertainties through the calculation chain using the stochastic sampling method. Output responses including assembly burnup, critical boron concentration, peak power, and cycle length were evaluated. The maximum uncertainty in assembly burnup is found to be 0.82% and the uncertainty of cycle length is found to be 36.2 days, both of which are mostly contributed by the nuclear data uncertainty.
It was shown in the current and previous studies that the uncertainty of certain core responses is not equally contributed by the uncertainty of three physics domains under consideration. However, this does not mean that the impact of input uncertainty should be propagated and quantified in an independent way; that is, the correlations of input variables due to the common source of uncertainties must be accounted for consistently. Taking the nuclear data uncertainty for example, it obviously introduces uncertainties to the predicted fuel composition along the burnup cycle, which in turn has impact on the neutronics parameters (i.e., fewgroup cross sections), the boundary condition of the thermalhydraulics calculation (i.e., heat flux), and the fuel modelling parameters (e.g., fuel thermal conductivity and gap conductance). All those variables are thus correlated and so are their uncertainties. Therefore, it is planned to develop a consistent uncertainty quantification method for the coupled core calculation by constructing the correlation among selected input parameters, which is represented by the components and data flow in dashed line in Figure 16.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors would like to acknowledge and recognise K. Zeng, N. Porter, and A. Toptan, current and former Ph.D. students, at the Reactor Dynamics and Fuel Modelling Group (RDFMG), Nuclear Engineering Department, North Carolina State University, who have contributed to the preparation of OECD/NEA LWR UAM benchmark specifications and development of methods as well as obtaining and analysing the results presented in this manuscript.
References
 K. Ivanov, M. Avramova, S. Kamerow et al., “Benchmark for uncertainty analysis in modelling (UAM) for the design, operation and safety analysis of LWRs,” in Specification and Support Data for Neutronics Cases (Phase I), vol. I, p. 7, NEA/NSC/DOC, Paris, France, 2013. View at: Google Scholar
 M. Avramova, “Benchmark for uncertainty analysis in modelling (UAM) for the design, operation and safety analysis of LWRs,” in Specification and Support Data for Neutronics Cases (Phase II), vol. II, NEA/NSC/DOC, Paris, France, 2017. View at: Google Scholar
 J. Hou, “Benchmark for uncertainty analysis in modelling (UAM) for the design, operation and safety analysis of LWRs,” in Specification and Support Data for Neutronics Cases (Phase III), vol. III, NEA/NSC/DOC, Paris, France, 2018. View at: Google Scholar
 M. Pusa, “Perturbationtheorybased sensitivity and uncertainty analysis with CASMO4,” Science and Technology of Nuclear Installations, vol. 2012, Article ID 157029, 11 pages, 2012. View at: Publisher Site  Google Scholar
 B. Rearden, “TSUNAMI1D: Control Module for OneDimensional CrossSection Sensitivity and Uncertainty”, Oak Ridge National Laboratory, Oak Ridge, TN, USA, 2011.
 B. Krzykacz, E. Hofer, and M. Kloos, “A software system for probabilistic uncertainty and sensitivity analysis of results from computer models,” in Proceedings of the International Conference Probabilistic Safety Assessment and Management (PSAMII), San Diego, CA, USA, January 1994. View at: Google Scholar
 M. Williams, “A statistical sampling method for uncertainty analysis with SCALE and XUSA,” Nuclear Technology, vol. 183, p. 515, 2012. View at: Google Scholar
 B. Rearden and M. Jessee, SCALE Code System, ORNL/TM2005/39, Version 6.2.1, UT Battelle, LLC, Oak Ridge Laboratory, Oak Ridge, TN, USA, 2016.
 W. Marshall, M. L. Williams, D. Wiarda et al., “Development and testing of neutron crosssection covariance data for SCALE 6.2,” in Proceedings of the International Conference on Nuclear Criticality Safety (ICNC 2015), Charlotte, NC, USA, September 2015. View at: Google Scholar
 A. Trkov, G. L. Molnár, Z. Révay et al., “Revisiting the ^{238}U thermal capture cross section and gammaray emission probabilities from ^{239}Np decay,” Nuclear Science and Engineering, vol. 150, no. 3, pp. 336–348, 2005. View at: Publisher Site  Google Scholar
 K. Shibata, T. Kawano, T. Nakagawa et al., “Japanese evaluated nuclear data library version 3 revision3: JENDL3.3,” Journal of Nuclear Science and Technology, vol. 39, no. 11, pp. 1125–1136, 2002. View at: Publisher Site  Google Scholar
 E. Castro, S. SánchezCervera, N. GarcíaHerranz, and D. Cuervo, “Impact of the homogenization level, nodal or pinbypin, on the uncertainty quantification with core simulators,” Progress in Nuclear Energy, vol. 104, pp. 218–228, 2018. View at: Publisher Site  Google Scholar
 I. Yankov, B. Collins, M. Klein et al., “A twostep approach to uncertainty quantification of core simulators,” Science and Technology of Nuclear Installations, vol. 2012, Article ID 767096, 9 pages, 2012. View at: Publisher Site  Google Scholar
 K. Zeng, J. Hou, K. Ivanov, and M. A. Jessee, “Uncertainty quantification and propagation of multiphysics simulation of the pressurized water reactor core,” Nuclear Technology, vol. 205, no. 12, pp. 1618–1637, 2019. View at: Publisher Site  Google Scholar
 K. Ivanov and M. Avramova, “Uncertainty analysis in modelling for Light water reactors consistent approach for multiscale modelling,” in Proceedings of the ANS Best Estimate Plus Uncertainty International Conference (BEPU 2018), Real Collegio, Lucca, Italy, May 2018. View at: Google Scholar
 G.K. Delipei, J. Garnier, J.C. Le Pallec, and B. Normand, “Uncertainty analysis methodology for multiphysics coupled rod ejection accident,” in Proceedings of the International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering (M&C), Portland, OR, USA, August 2019. View at: Google Scholar
 K. Zeng, J. Hou, K. Ivanov, and J. Matthew, “Uncertainty analysis of Light water reactor core simulations using statistic sampling method,” in Proceedings of the International Conference on Mathematics & Computational 2017, Jeju, Korea, April 2017. View at: Google Scholar
 N. W. Porter, M. N. Avramova, and V. A. Mousseau, “Uncertainty quantification study of CTF for the OECD/NEA LWR uncertainty analysis in modeling benchmark,” Nuclear Science and Engineering, vol. 190, no. 3, pp. 271–286, 2018. View at: Publisher Site  Google Scholar
 M. Avramova, “Developments in thermalhydraulic subchannel modeling for whole core multiphysics simulations,” Nuclear Engineering and Design, vol. 358, p. 110387, 2020. View at: Publisher Site  Google Scholar
 T. Downar, Y. Xu, and V. Seker, “PARCS–U.S. NRC core neutronics simulator, UMNERS090001,” 2013. View at: Google Scholar
 A. Wysocki, A. Ward, A. Manera et al., “The modeling of advanced BWR fuel designs with the NRC fuel depletion codes PARCS/PATHS,” Nuclear Technology, vol. 190, no. 3, pp. 323–335, 2015. View at: Publisher Site  Google Scholar
 B. Adams, Dakota, A Multilevel Parallel ObjectOriented Framework for Design Optimization, Parameter Estimation, Uncertainty Quantification, and Sensitivity Analysis: Version 6.5 User’s Manual, Sandia National Laboratory, Albuquerque, NM, USA, 2016.
 S. S. Wilks, “Determination of sample sizes for setting tolerance limits,” The Annals of Mathematical Statistics, vol. 12, no. 1, pp. 91–96, 1941. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2020 Jason Hou et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.