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 Best-Estimate 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 stand-alone 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. High-to-Low (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 time-dependent 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 higher-fidelity models implemented in fuel performance codes and the rather simplified models implemented in thermal-hydraulics codes, to be used for coupling with neutronics in Phase III is presented. Similarly, the uncertainty quantification on thermal-hydraulic 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 up-to-date 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 best-estimate 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 best-estimate 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 best-estimate 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, thermal-hydraulics, 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 TMI-1 reactor related test cases. Section 3 summarizes one of the concluding activities of the stand-alone 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 thermal-hydraulics 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 first-order uncertainty propagation formula or the so-called “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 chi-square distribution with N − 1 degree of freedom:where is the sample variance. For a two-sided 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 two-sided 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. TMI-1 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 (TMI-1) reactor, which is used as an example to illustrate the developments and implementations performed with the LWR UAM benchmark framework.
The TMI-1 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 TMI-1 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 (GdO2 + UO2), 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 stand-alone 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 full-core performance parameters). The sources of the input uncertainties considered in Phase I exercises include the neutron cross section data, supplemented by the VCMs, and as-built 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 44-group library distributed with SCALE 6.0 and SCALE 6.1 and the updated library available in SCALE 6.2 which is based on 56-group and 252-group structures [8]. The former contains uncertainty data for 401 materials with important isotopes taken from high-fidelity nuclear data evaluations including ENDF/B-VII.0, ENDF/B-VI, and JENDL-3.3. The latter is based on the EDNF/B-VII.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.
3.1. Cell Physics
For each test case in Exercise I-1, the following results are requested from the participants: calculated k-inf and its associated uncertainty, the top five neutron-nuclide reactions that contribute the most uncertainty to k-inf, and covariances of selected one-group cross sections generated in the pin cell calculation.
A two-dimensional (2D) model adopted from the TMI-1 reactor was chosen as the representative PWR test problem. It is fuelled with 4.85% enriched UO2 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 k-inf and associated uncertainties. The nominal value of k-inf 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 k-inf 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 k-inf and the results indicate that the choice of the covariance library strongly impacts the RSD of k-inf, as shown in Figure 4, while other parameters have limited influence. The average RSD of k-inf 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 JENDL-4.0 covariance libraries and the RSDs are ∼0.5%. The ENDF/B-VII.1 covariance data yields the highest k-inf 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 k-inf 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 nu-bar (average number of neutrons per fission reaction ) of 235U and 239Pu, with the first one relevant to the case study shown in this paper. The uncertainty of 235U nu-bar in the thermal range increases from 0.31% in SCALE 6.1 to 0.7% in ENDF/B-VII.1, which is responsible for differences exhibited in the covariance testing in the low-enriched, 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 JENDL-3.3 [8]. This explains why the calculated RSD of k-inf using SCALE 6.2 data is slightly higher than that using SCALE 6.0/SCALE 6.1, while the value corresponding to ENDF/B-VII.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 nuclide-reaction pair to the predicted k-inf uncertainties by sorting them from greatest to least variance fraction. For the TMI-1 HZP case, 13 sets of submitted results include such information and Figure 5 shows the occurrence of various nuclide-reaction pairs as the top 5 contributors. Although a certain degree of diversity can be found in the ranking as it includes up to 10 nuclide-reaction pairs, some reactions dominate the contribution to the uncertainty of k-inf, such as 238U capture, 235U nu-bar, and 235U 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 k-inf is quite sensitive to the 238U capture cross section, especially in the unresolved resonance regions, while the evaluated cross sections exhibit large uncertainties [10]. This is why, on one hand, 238U 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, 235U nu-bar tops the ranking if the SCALE 6.2 or ENDF/B-VII.1 covariance libraries are used, which again confirms the evolution of the 235U nu-bar uncertainty as described above. The 235U 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 nuclide-reaction pairs is depicted in Figure 6. It is worth mentioning that some nuclide-reaction pairs only associate with the ENDF/B-VII.1 covariance library, including the 1H capture, 238U nu-bar, and 238U elastic scattering, which can also be explained by the difference in covariance libraries. For example, the 1H thermal capture in SCALE covariance library is adopted from JENDL 3.3 [11], which is lower than that in ENDF/B-VII.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 few-group 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 k-inf 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 k-inf 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 235U nu-bar.
(a)
(b)
The uncertainties of the homogenized two-group nu-fission cross sections are compared based on the VCM, as shown in Figure 8. Again, the higher nuclear data uncertainty of the nu-bar 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 nu-fission 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 235U nu-bar 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 k-inf with respect to the 235U nu-bar and fission cross section [12]. The comparison reflects the fact that the nominal value of nu-fission 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 k-inf, flux, and group constants (cross sections, diffusion coefficients, and assembly discontinuity factors) in the two-group 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 k-inf and Group 2 nu-fission cross section. On the other hand, some unique features are also observed only in one of the cases. Case 12 shows strong cross-correlations 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 k-inf 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 few-group cross sections evaluated in Exercise I-2 are to be propagated through the stand-alone neutronics core calculation in Exercise I-3 to the parameters of interest such as the core k-eff 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 few-group homogenized constants generated in the lattice calculation is used to evaluate the uncertainty of core responses. The one-step approach relies on the stochastic sampling method on both lattice and core levels and is named as such because it involves a one-to-one connection between lattice calculations to generate the few-group cross section library and the core calculation that read this library as input [13]. The two-step 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 TMI-1 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 one-step and two-step approaches mentioned above, such as in [14], except for Cases 11 and 35, where the direct full-core 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 full-core 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 Thermal-Hydraulics 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 time-dependent 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 II-2 considers time-dependent 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 II-1 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 thermal-hydraulics codes, to be used for coupling with neutronics tools in Phase III [15].
High-fidelity fuel performance codes are used to model a single fuel rod, using detailed technological and operation data (geometry, enrichment, burnup, etc.) and high-fidelity models for fission gas release, cladding corrosion, swelling, and so forth. These codes are used in Exercise II-1, while low-fidelity models, embedded in thermal-hydraulics codes, are usually preferred for full-core computation in Phase III. In this way Exercise II-1 propagates uncertainties in fuel physics consistently using high-fidelity 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 high-fidelity fuel performance codes for the PWR depletion test case based on TMI-1 fuel rod design. Exercise II-1 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 thermal-hydraulics codes to predict fuel temperature as compared to fuel performance codes.
The so-called Doppler temperature is the main fuel physical feedback quantity used in full-core 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 thermal-hydraulics codes) or for an individual rod in a pin-by-pin (or cell-by-cell) calculation (in fuel performance codes). In both cases, the fuel temperature is a result of a one-dimension 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 thermal-hydraulic conditions as solved by the core thermal-hydraulics model. The rod-to-coolant heat transfer coefficient is also given by thermal-hydraulics 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 high-fidelity 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 II-1 with higher-fidelity 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 high-fidelity fuel performance codes—(i) the mean value and (ii) associated uncertainty—both as functions of the following parameters: for cladding conductivity kc(T), as function of temperature (T); for fuel conductivity kf(Bu, T), as function of burnup (Bu) and temperature (T); and for gap conductance, hg(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 high-fidelity fuel models for lumped simplified lower-fidelity fuel models. The lookup table representation concept [2] is appropriate for steady state and cycle depletion calculations while for transient/accident simulations the high-to-low 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. Thermal-Hydraulics
As for fuel thermal and mechanical behaviour, there might be different thermal-hydraulic models used between Phases II and III. The uncertainty quantification on thermal-hydraulic models is established on a relatively small scale, that is, rod bundles, in Exercise II-3, while these results will be used in Phase III at the core scale. Detailed two-phase flow models, either at subchannel level or at Computational Fluid Dynamics (CFD) level, are used to model a single fuel bundle in Exercise II-3, while rather simplified models at fuel assembly level are utilized for full-core thermal-hydraulics in Phase III [15].
Considering that the feedback effect of thermal-hydraulics 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: single-phase turbulent flow, nucleate boiling, or film boiling (post-CHF 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 two-phase 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., four-equation model versus two-fluid 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 TMI-1 PWR assembly at normal operation conditions as part of Exercise II-3, 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 TMI-1 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 single-phase 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 single-phase flow is simpler and better understood than two-phase modelling, which leads to smaller overall uncertainty in thermal-hydraulic 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 in-core structure modelling [19]. These CFD-informed models are developed using a physics-based approach which utilizes high-fidelity data to inform and calibrate a lower-fidelity simulation. It is important to consider only once input uncertainties. For example, the spacer grid effects (position and mixing) are investigated in Exercise II-3 for a given type of fuel bundle—in this case, TMI-1 PWR bundle. The resulting uncertainty will be considered later in Phase III in the overall uncertainty on investigated output parameters. The methodology utilizes CFD-based data for: grid-directed cross-flow model; grid-enhanced turbulent mixing model; grid-enhanced heat transfer model; and spacer grid pressure losses model. This data while utilized in lower-fidelity core thermal-hydraulics modelling as part of core multiphysics calculations in Phase III is supplemented with uncertainties propagated through the high-fidelity 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, thermal-hydraulics, 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 III-1, 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 III-1. 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 thermal-hydraulics.
The neutronics uncertainties are embedded in the two-group 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 kf, gap conductance hg, and cladding thermal conductivity kc. In general, the input uncertainties for the core thermal-hydraulic 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 TMI-1 equilibrium cycle depletion case was specified in the LWR UAM benchmark, Exercise III-1, 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 TMI-1 core was modelled using PARCS [20] with the fixed mesh size. The core is radially partitioned into 177 nodes based on the one-node-per-assembly 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 two-group homogenized cross sections were generated based on the ENDF/B-VII.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 Thermal-Hydraulics Solver) code [21] to provide thermal-hydraulics feedback for steady state and depletion calculation via the tight coupling with PARCS. The axial nodalization and spatial mapping (overlays) of thermal-hydraulics (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 56-group covariance library and generate random samples of nuclear data for the subsequent lattice calculation. The few-group constants were then generated and converted by the lattice physics code Polaris [8]. On the thermal-hydraulics 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 few-group constants with one sample of fuel modelling and thermal-hydraulics 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 python-based 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 95th 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 thermal-hydraulics 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. Thermal-hydraulic 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 thermal-hydraulics 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 Best-Estimate 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 TMI-1 reactor. The stand-alone 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 k-inf (k-eff) uncertainty strongly depends on the VCM and the results generated using the VCM of ENDF/B-VII.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 235U nu-bar uncertainty in the thermal energy range is responsible for the disagreement. As a follow-up to the current work, similar analyses will be applied to benchmark cases belonging to other reactor types, including BWR, VVER, and Gen-III 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 II-2 to parameterize the neutronics uncertainties in a form of few-group libraries to be propagated from Phase II to Phase III. Exercises II-1 and II-3 are used to characterize the uncertainties of high-fidelity fuel and thermal-hydraulics 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 thermal-hydraulics. For the three selected fuel physics quantities, high-to-low-fidelity 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 thermal-hydraulic quantities where common cross-cutting 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 full-core calculation. For such propagation, the experience accumulated in Exercise II-3 test problems is utilized. Finally, the multiphysics exercises of Phase III will result in safety and design quantities with the best-estimate 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., few-group cross sections), the boundary condition of the thermal-hydraulics 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.