# Lognormal-Based Sampling for Fission Product Yields Uncertainty Propagation in Pebble-Bed HTGR

Authors: Yizhen Wang; Menglei Cui; Jiong Guo; Jinlin Niu; Yingjie Wu; Baokun Liu; Fu Li
Journal: Science and Technology of Nuclear Installations
 (2020)

---

## Abstract

Uncertainty analyses of fission product yields are indispensable in evaluating reactor burnup and decay heat calculation credibility. Compared with neutron cross section, there are fewer uncertainty analyses conducted and it has been a controversial topic by lack of properly estimated covariance matrix as well as adequate sampling method. Specifically, the conventional normal-based sampling method in sampling large uncertainty independent fission yields could inevitably generate nonphysical negative samples. Cutting off these samples would introduce bias into uncertainty results. Here, we evaluate thermal neutron-induced U-235 independent fission yields covariance matrix by the Bayesian updating method, and then we use lognormal-based sampling method to overcome the negative fission yields samples issue. Fission yields uncertainty contribution to effective multiplication factor and several fission products’ atomic densities at equilibrium core of pebble-bed HTGR are quantified and investigated. The results show that the lognormal-based sampling method could prevent generating negative yields samples and maintain the skewness of fission yields distribution. Compared with the zero cut-off normal-based sampling method, the lognormal-based sampling method evaluates the uncertainty of effective multiplication factor and atomic densities are larger. This implies that zero cut-off effect in the normal-based sampling method would underestimate the fission yields uncertainty contribution. Therefore, adopting the lognormal-based sampling method is crucial for providing reliable uncertainty analysis results in fission product yields uncertainty analysis.


---

## Body

## 1. Introduction

Reactor design and safety analysis rely on accurate calculations of system responses with properly evaluated uncertainties. There has been an increasing need for evaluating the credibility of reactor safety. Pebble-bed high temperature gas-cooled reactor (pebble-bed HTGR) is a multiphysics nonlinear coupled system, including neutron transport and complex heat transfer hydraulics behaviour [1]. In order to systematically and thoroughly investigate the uncertainties propagation in pebble-bed HTGR, an IAEA Coordinated Research Plan (CRP) [2, 3] has been initiated after the start of OECD/NEA UAM-LWR [4]. Recent advances regarding uncertainty propagation analysis in pebble-bed HTGR mainly concern the nuclear cross section uncertainties propagation in reactor neutronic calculations [5–9]. As pebble-bed HTGR allows fuels recirculation during fuel cycles and adopts higher fuel enrichment (8.5 wt.%), fuels usually could achieve larger burnup values, and then fission product yields could be nonnegligible uncertainty sources in reactor burnup and decay heat calculations. Their uncertainties contributions to important reactor burnup responses need to be considered properly for evaluating the credibility of rector safety-related quantities of interest (QoI), e.g., maximum fuel pebble temperature.Fission product yield describes the fraction of a certain fission product produced per fission. During the measurements of fission product yields, correlated errors or covariances may exist when using the same equipment or methods [10]. However, they are ignored in evaluated nuclear data library. Also, self-consistent fission yields data set should follow several physical constraints such as binary fission, mass conservation, and charge conservation [11]. These constraints could introduce covariances between fission yields data. As the fission yields’ covariances in current releases of evaluated nuclear data libraries are still absent, e.g., ENDF/B-VII.1, many methods are developed to estimate these covariances information based on the imposed physical constraints. The Bayesian updating method is widely used in data assimilation, data adjustment, and model fitting problems. It refines parameters by taking both the prior information about those parameters and the likelihood which refers to new data into consideration [10, 12]. It allows estimating the covariance matrix of fission yields by sequentially introducing the above physical constraints. In the domain of fission yields adjustment, it is introduced by Kawano and Chadwick [13] to update Pu-239 fission yields with chain yields to reduce the independent fission yields discrepancy in ENDF/B-VII.1. Pigni et al. [14] expand it to involve cumulative fission yields into covariances estimation. The difference between chain yields-based updating and cumulative fission yields-based updating is further investigated by Fiorito et al. [15, 16]. Based on the provided independent and cumulative fission yields uncertainties information in ENDF/B-VII.1, this work adopts the Bayesian updating method to estimate the independent fission yields covariances.Sampling-based methods for uncertainty analysis or stochastic UQ methods [17] require properly perturbed samples to provide reliable uncertainty analysis results of QoI. As it is observed from the evaluated nuclear data library, independent fission yields generally have larger uncertainties. Random sampling on these yields under normal distribution could generate nonphysical negative samples. Cutting off these negative yield samples by setting them to zero could artificially omit some information from the original fission yields distribution, resulting in biased uncertainty analysis results. This zero cut-off effect on quantified uncertainty has not been well studied. This raises question whether normal distribution is sufficient for describing inherently positive random variables with large uncertainties only given their mean values and covariance matrix. Smith et al. [18] propose to replace normal distribution with lognormal distribution by the principle of maximum entropy [10], and Žerovnik et al. [19, 20] investigate this method in the sampling resonance parameters where negative samples problem was encountered as in fission yields. This work proposes an implementation of the lognormal-based sampling method in fission product yields sampling.The present work is organized as follows: Section2 describes the nomenclature of fission product yields and the burnup calculation of pebble-bed HTGR. An implemented stochastic UQ method for fission yields uncertainty propagation is described in Section 2.3. The Bayesian updating method and the lognormal-based sampling method are detailed in Section 3. Finally, results of fission yields uncertainty contributions to effective multiplication factor and several important fission products atomic densities are provided and discussed in Section 4.


## 2. Model Description and Uncertainty Propagation

### 2.1. ENDF/B-VII.1 Fission Product Yields Sublibrary

Fission product yield characterizes the fraction of a particular fission product nuclide produced per fission. A compound nucleus is formed when a fissile nucleus is bombarded by an incident neutron. As its energy overcomes the fission barrier, this compound nucleus could undergo fission. A brief description of the fission process is illustrated (see Figure1). After the scission of compound nucleus, primary fission fragments are produced and they would undergo deexcitation by releasing prompt neutrons due to their high neutron to proton ratios. After the emission of prompt neutrons, the remaining fission fragments referred to as fission products would undergo β decay, isomeric transition, or particle emission along their corresponding decay chain and finally reach stable nuclides. Each fission product is identified by its mass number A, charge number Z, and isomeric state I and is denoted as the triplet A,Z,IFigure 1
Neutron induced fission process [21]. The fission products refer to the fission fragments after the emission of prompt neutrons. Independent fission yields (IFYs) characterize the fraction of a fission product produced before any radioactive decay, whereas cumulative fission yields (CFYs) describe the fraction of that produced product over all time after a fission.A detailed description about the nomenclature of fission product yields could be found in [11] and they are briefly summarized as follows. IFYs and CFYs determine the fraction of a fission product at different stages in the fission process. IFYs denoted as yA,Z,I are the fraction of a fission product produced directly from one fission after the emission of prompt neutrons but prior to any radioactive decays. Because IFYs are produced before any radioactive decay in the fission system, they should be subject to the physical constraints of fission system, e.g., binary fission, conservation of mass, and charge number. CFYs denoted as cA,Z,I determine the total fraction of a fission product produced over all time after one fission. It involves not only the direct production from fission but also the contributions from the decay of other products.The current releases of ENDF/B-VII.1 fission yield sublibrary provide fission yields data for 31 fission actinides from Th-227 to Fm-255. Though energy-dependence issues within fission spectrum are highlighted in current releases of evaluated nuclear data library and neutron induced Pu-239 fission yields at 2.0 MeV are supplemented to allow users to linearly interpolate yields between 0.5 MeV and 2.0 MeV for high accuracy purpose [22], other fission actinide fission yields data are taken directly from ENDF/B-VI evaluated by England and Rider [23] in 1993. Three fission systems for U-235 are evaluated with respect to incident neutron energy, namely, 0.0253 eV thermal energy, 0.5 MeV fission spectrum energy, and 14.0 MeV high energy. IFYs and CFYs are evaluated for 1,247 fission products in thermal neutron induced U-235 fission yield (see Figure 2). The relationship between IFYs and CFYs [11] is referred to as (1), where bA′,Z′,I′⟶A,Z,I is the branching ratio:(1)cA,Z,I=yA,Z,I+∑A′,Z′,I′bA′,Z′,I′⟶A,Z,IcA′,Z′,I′.Figure 2
Thermal neutron induced U-235 fission yield data in ENDF/B-VII.1. The natural logarithm values of fission yields are plotted. (a) Independent yields. (b) Cumulative yields.
(a)(b)It could be found that most IFYs appear in the upper region ofβ-stability line and they are most likely to undergo β− decay to reach a stable state. As CFYs involve the production of a certain fission product from the decay of other fission products as shown in (1), the peaks of CFYs distribution in neutron-charge number figure tend to be closer to the β-stability line (see Figure 2).The evaluation of fission yields data requires a combined work of experimental measurements and theoretical model predictions. It is natural for the evaluated fission yields possessing uncertainties originated from measurement errors and theoretical model parameters uncertainties. Although England and Rider provide the uncertainties (standard deviation) of each fission product yield in their original work, covariances information between fission yields has not been provided since then. Those covariances information is crucial for representing the physical constraints imposed on IFYs, and they should be estimated properly in order to generate self-consistent IFYs perturbations. This work focuses on the propagation of thermal neutron induced U-235 fission yields uncertainties to burnup simulation of pebble-bed HTGR based on ENDF/B-VII.1. The estimation of covariances information will be detailed in Section3.1.


### 2.2. Pebble-Bed HTGR Burnup Model and Built-In Fission Yields Analysis

Pebble-bed HTGR core (see Figure3(a)) consists of spherical fuel elements or fuel pebbles. Each of these pebbles is composed of a spherical graphite matrix in the centre where thousands of small coated particles (known as TRISO particles) are embedded. These particles contain UO2 kernel in the centre with four structural coating layers surrounding it (see Figure 3(b)). During reactor operation, these fuel pebbles are consistently flowing downward from the top of the core to the bottom and are irradiated at different core spectrum regions randomly. Fuel recirculation is a characterized fuel cycling procedure adopted in pebble-bed HTGR which is different with that applied in Light Water Reactor (LWR). Such recirculation allows fresh fuel pebbles being loaded into the core and spent fuel pebbles being discharged online without shutting reactor down. More importantly, this recirculation permits fuel pebbles running through core multiple times before they are finally being discharged. Because of the fuels recirculation, there exist running-in phase and equilibrium core states. The equilibrium core state refers to the nuclei compositions inside the core kept unchanged with time and therefore effective multiplication factor being stable at a certain value. This could give a more flattened power distribution across the core and higher average discharge burnup value. The V.S.O.P. computer code system [26] is developed to perform burnup calculation of pebble-bed HTGR by simulating the fuels recirculation process stepwise and conduct spectrum calculation online at each spectrum region inside the core. A detailed description of this simulation process could be found in these articles [27].Figure 3
Pebble-bed HTGR core. (a) Core geometrical of PBR250 design [24]. (b) Fuel pebbles [25].
(a)(b)The built-in fission product chain in V.S.O.P. code involves 44 fission products and among these 44 fission products’ data, 14 are taken as IFYs, while 30 are taken as CFYs. These data are taken from ENDF/B-IV and ENDF/B-V. An additional “nonsaturating” fission product is evaluated to account for the sum of many lumped fission yields which are not explicitly included in the chain [26]. The comparison between these built-in fission yields and replacing them with the current releases in ENDF/B-VII.1 is conducted to examine the availability of V.S.O.P. code for fission yields uncertainty propagation. This investigation is conducted on HTR-PM [28] with 8.7% fuel enrichment (while 8.5% enrichment is applied in actual design).The impact of each built-in fission yield onkeff at equilibrium core state is investigated individually by replacing them with ENDF/B-VII.1. It should be noted that built-in fission yields library in V.S.O.P. includes a combination of IFYs and CFYs, and they are presented separately (see Tables 1 and 2).Table 1
Comparison between V.S.O.P. built-in yields and ENDF/B-VII.1 inkeff prediction (IFYs).
IndexFission productFission yieldskeffV.S.O.P.ENDF/B-VII.1Relative difference to ENDF/B-VII.1 (%)V.S.O.P.ENDF/B-VII.1Difference to ENDF/B-VII.1 (pcm)1Rh-1031.8580E − 116.3796E − 1396.571.010271.0102702Pd-1059.8300E − 130.0000E + 00100.001.010271.0102703Xe-1311.5400E − 081.4199E − 0990.781.010271.0102704Cs-1335.0800E − 077.9194E − 0998.441.010271.0102705Cs-1343.5700E − 073.8547E − 0889.201.010271.0102706Nd-1439.5000E − 134.7997E − 1494.951.010271.0102707Pm-148m7.4900E − 098.0994E − 1198.921.010271.0102708Pm-148g5.7300E − 084.4497E − 1199.921.010271.0102709Sm-1470.0000E + 000.0000E + 000.001.010271.01027010Sm-1486.9500E − 131.6399E − 1497.641.010271.01027011Sm-1490.0000E + 001.7099E − 12—1.010271.01027012Sm-1510.0000E + 004.7497E − 09—1.010271.01027013Eu-1541.6300E − 089.6993E − 1094.051.010271.01027014Gd-1554.4100E − 114.0797E − 1290.751.010271.010270Table 2
Comparison between V.S.O.P. built-in yields and ENDF/B-VII.1 inkeff prediction (CFYs).
IndexFission productFission yieldskeffV.S.O.P.ENDF/B-VII.1Relative difference to ENDF/B-VII.1 (%)V.S.O.P.ENDF/B-VII.1Difference to ENDF/B-VII.1 (pcm)1Xe-1356.6023E − 026.5385E − 020.971.010271.01044−172FP-449.4760E − 019.4760E − 010.001.010271.0102703Xe-1366.2701E − 026.3127E − 02−0.681.010271.0102704Kr-835.3076E − 035.3620E − 03−1.021.010271.0102705Zr-956.4678E − 026.5027E − 02−0.541.010271.0102706Mo-951.6410E − 066.5029E − 02−3.96×10101.010271.009111167Mo-975.9600E − 025.9968E − 02−0.621.010271.0102708Tc-996.1284E − 026.1087E − 020.321.010271.01028−19Ru-1015.0501E − 025.1725E − 02−2.421.010271.01026110Ru-1033.1411E − 023.0309E − 023.511.010271.01042−1511Rh-1051.0199E − 029.6416E − 035.471.010271.01030−312Pd-1087.1032E − 045.4125E − 0423.801.010271.01028−113Ag-1092.9903E − 043.1221E − 04−4.411.010271.01027014Cd-1131.2425E − 041.4038E − 04−12.981.010271.01027015I-1312.8325E − 022.8907E − 02−2.051.010271.01022516Xe-1336.7859E − 026.6991E − 021.281.010271.01032−517Pr-1415.8929E − 025.8470E − 020.781.010271.01028−118Pr-1435.9710E − 025.9558E − 020.251.010271.01029−219Nd-1445.4523E − 025.4996E − 02−0.871.010271.01027020Nd-1453.9339E − 023.9334E − 020.011.010271.01027021Nd-1462.9912E − 022.9969E − 02−0.191.010271.01027022Pm-1472.2701E − 022.2467E − 021.031.010271.01035−823Pm-1491.0888E − 021.0816E − 020.591.010271.01031−424Sm-1505.4130E − 062.9998E − 0794.461.010271.01027025Pm-1514.2044E − 034.1877E − 030.401.010271.01028−126Sm-1522.7057E − 032.6691E − 031.351.010271.01029−227Eu-1531.6264E − 031.5828E − 032.681.010271.01029−228Eu-1553.3025E − 043.2136E − 042.691.010271.01028−129Gd-1561.3517E − 041.4853E − 04−9.881.010271.01027030Gd-1576.4651E − 056.1506E − 054.861.010271.010270IFYs are evaluated by subtracting the total contributions of its precursors from experimental measured CFYs. With the improvement of CFYs measurements, the evaluated IFYs become more precise. It could be seen from the table that IFYs in ENDF/B-VII.1 are lower than the built-in fission yields used in V.S.O.P. Except the large discrepancy in the fission yield of Mo-95, all the impacts from replacing fission yields are lower than 20 pcm. The overall impact is 67 pcm (see Table3) when all the yields are replaced without FP-44. The difference is acceptable in effective multiplication factor calculations when substituting built-in V.S.O.P. fission yields with ENDF/B-VII.1 fission yields. The V.S.O.P. burnup model is further used to conduct fission product yields uncertainty propagation as described in Section 2.3.Table 3
Reference calculation between built-in fission yields of V.S.O.P. and ENDF/B-VII.1 (all substitution without FP-44).
keffDifference to ENDF/B-VII.1 (pcm)Built-in V.S.O.P. fission yieldsENDF/B-VII.1 fission product yields (without update)1.010271.0096067


### 2.3. Uncertainty Quantification Scheme

The HTR-PM [28] reactor core is modelled in V.S.O.P. computer code system to analyse the uncertainty propagation of fission yields in this work. 15 times recirculation of fuel is adopted and the average discharge burnup value is around 90,210MW⋅d/tU with fresh fuel having 8.5 wt.% enrichment. As fuel recirculation tightly couples the neutronics and burnup calculation spatially inside the core, it is difficult to separate the uncertainty propagation step by step. Stochastic UQ method is used to investigate the uncertainty propagation in equilibrium core state. An uncertainty propagation scheme is proposed in this work (see Figure 4).Figure 4
Flow chart of fission product yields uncertainty propagation.Two sampling methods are implemented in this work, namely, normal-based sampling and lognormal-based sampling. Different from normal-based sampling, lognormal-based sampling requires a lognormal transformation of the original mean vector and covariance matrix. When the IFYs samples are generated, their corresponding CFYs are calculated and combine them to form self-consistent fission yield samples. These combined IFYs and CFYs samples are propagated to V.S.O.P. HTR-PM model for further uncertainty analysis. Detailed Bayesian updating method description and lognormal-based sampling procedures will be introduced in Section3.


## 2.1. ENDF/B-VII.1 Fission Product Yields Sublibrary

Fission product yield characterizes the fraction of a particular fission product nuclide produced per fission. A compound nucleus is formed when a fissile nucleus is bombarded by an incident neutron. As its energy overcomes the fission barrier, this compound nucleus could undergo fission. A brief description of the fission process is illustrated (see Figure1). After the scission of compound nucleus, primary fission fragments are produced and they would undergo deexcitation by releasing prompt neutrons due to their high neutron to proton ratios. After the emission of prompt neutrons, the remaining fission fragments referred to as fission products would undergo β decay, isomeric transition, or particle emission along their corresponding decay chain and finally reach stable nuclides. Each fission product is identified by its mass number A, charge number Z, and isomeric state I and is denoted as the triplet A,Z,IFigure 1
Neutron induced fission process [21]. The fission products refer to the fission fragments after the emission of prompt neutrons. Independent fission yields (IFYs) characterize the fraction of a fission product produced before any radioactive decay, whereas cumulative fission yields (CFYs) describe the fraction of that produced product over all time after a fission.A detailed description about the nomenclature of fission product yields could be found in [11] and they are briefly summarized as follows. IFYs and CFYs determine the fraction of a fission product at different stages in the fission process. IFYs denoted as yA,Z,I are the fraction of a fission product produced directly from one fission after the emission of prompt neutrons but prior to any radioactive decays. Because IFYs are produced before any radioactive decay in the fission system, they should be subject to the physical constraints of fission system, e.g., binary fission, conservation of mass, and charge number. CFYs denoted as cA,Z,I determine the total fraction of a fission product produced over all time after one fission. It involves not only the direct production from fission but also the contributions from the decay of other products.The current releases of ENDF/B-VII.1 fission yield sublibrary provide fission yields data for 31 fission actinides from Th-227 to Fm-255. Though energy-dependence issues within fission spectrum are highlighted in current releases of evaluated nuclear data library and neutron induced Pu-239 fission yields at 2.0 MeV are supplemented to allow users to linearly interpolate yields between 0.5 MeV and 2.0 MeV for high accuracy purpose [22], other fission actinide fission yields data are taken directly from ENDF/B-VI evaluated by England and Rider [23] in 1993. Three fission systems for U-235 are evaluated with respect to incident neutron energy, namely, 0.0253 eV thermal energy, 0.5 MeV fission spectrum energy, and 14.0 MeV high energy. IFYs and CFYs are evaluated for 1,247 fission products in thermal neutron induced U-235 fission yield (see Figure 2). The relationship between IFYs and CFYs [11] is referred to as (1), where bA′,Z′,I′⟶A,Z,I is the branching ratio:(1)cA,Z,I=yA,Z,I+∑A′,Z′,I′bA′,Z′,I′⟶A,Z,IcA′,Z′,I′.Figure 2
Thermal neutron induced U-235 fission yield data in ENDF/B-VII.1. The natural logarithm values of fission yields are plotted. (a) Independent yields. (b) Cumulative yields.
(a)(b)It could be found that most IFYs appear in the upper region ofβ-stability line and they are most likely to undergo β− decay to reach a stable state. As CFYs involve the production of a certain fission product from the decay of other fission products as shown in (1), the peaks of CFYs distribution in neutron-charge number figure tend to be closer to the β-stability line (see Figure 2).The evaluation of fission yields data requires a combined work of experimental measurements and theoretical model predictions. It is natural for the evaluated fission yields possessing uncertainties originated from measurement errors and theoretical model parameters uncertainties. Although England and Rider provide the uncertainties (standard deviation) of each fission product yield in their original work, covariances information between fission yields has not been provided since then. Those covariances information is crucial for representing the physical constraints imposed on IFYs, and they should be estimated properly in order to generate self-consistent IFYs perturbations. This work focuses on the propagation of thermal neutron induced U-235 fission yields uncertainties to burnup simulation of pebble-bed HTGR based on ENDF/B-VII.1. The estimation of covariances information will be detailed in Section3.1.


## 2.2. Pebble-Bed HTGR Burnup Model and Built-In Fission Yields Analysis

Pebble-bed HTGR core (see Figure3(a)) consists of spherical fuel elements or fuel pebbles. Each of these pebbles is composed of a spherical graphite matrix in the centre where thousands of small coated particles (known as TRISO particles) are embedded. These particles contain UO2 kernel in the centre with four structural coating layers surrounding it (see Figure 3(b)). During reactor operation, these fuel pebbles are consistently flowing downward from the top of the core to the bottom and are irradiated at different core spectrum regions randomly. Fuel recirculation is a characterized fuel cycling procedure adopted in pebble-bed HTGR which is different with that applied in Light Water Reactor (LWR). Such recirculation allows fresh fuel pebbles being loaded into the core and spent fuel pebbles being discharged online without shutting reactor down. More importantly, this recirculation permits fuel pebbles running through core multiple times before they are finally being discharged. Because of the fuels recirculation, there exist running-in phase and equilibrium core states. The equilibrium core state refers to the nuclei compositions inside the core kept unchanged with time and therefore effective multiplication factor being stable at a certain value. This could give a more flattened power distribution across the core and higher average discharge burnup value. The V.S.O.P. computer code system [26] is developed to perform burnup calculation of pebble-bed HTGR by simulating the fuels recirculation process stepwise and conduct spectrum calculation online at each spectrum region inside the core. A detailed description of this simulation process could be found in these articles [27].Figure 3
Pebble-bed HTGR core. (a) Core geometrical of PBR250 design [24]. (b) Fuel pebbles [25].
(a)(b)The built-in fission product chain in V.S.O.P. code involves 44 fission products and among these 44 fission products’ data, 14 are taken as IFYs, while 30 are taken as CFYs. These data are taken from ENDF/B-IV and ENDF/B-V. An additional “nonsaturating” fission product is evaluated to account for the sum of many lumped fission yields which are not explicitly included in the chain [26]. The comparison between these built-in fission yields and replacing them with the current releases in ENDF/B-VII.1 is conducted to examine the availability of V.S.O.P. code for fission yields uncertainty propagation. This investigation is conducted on HTR-PM [28] with 8.7% fuel enrichment (while 8.5% enrichment is applied in actual design).The impact of each built-in fission yield onkeff at equilibrium core state is investigated individually by replacing them with ENDF/B-VII.1. It should be noted that built-in fission yields library in V.S.O.P. includes a combination of IFYs and CFYs, and they are presented separately (see Tables 1 and 2).Table 1
Comparison between V.S.O.P. built-in yields and ENDF/B-VII.1 inkeff prediction (IFYs).
IndexFission productFission yieldskeffV.S.O.P.ENDF/B-VII.1Relative difference to ENDF/B-VII.1 (%)V.S.O.P.ENDF/B-VII.1Difference to ENDF/B-VII.1 (pcm)1Rh-1031.8580E − 116.3796E − 1396.571.010271.0102702Pd-1059.8300E − 130.0000E + 00100.001.010271.0102703Xe-1311.5400E − 081.4199E − 0990.781.010271.0102704Cs-1335.0800E − 077.9194E − 0998.441.010271.0102705Cs-1343.5700E − 073.8547E − 0889.201.010271.0102706Nd-1439.5000E − 134.7997E − 1494.951.010271.0102707Pm-148m7.4900E − 098.0994E − 1198.921.010271.0102708Pm-148g5.7300E − 084.4497E − 1199.921.010271.0102709Sm-1470.0000E + 000.0000E + 000.001.010271.01027010Sm-1486.9500E − 131.6399E − 1497.641.010271.01027011Sm-1490.0000E + 001.7099E − 12—1.010271.01027012Sm-1510.0000E + 004.7497E − 09—1.010271.01027013Eu-1541.6300E − 089.6993E − 1094.051.010271.01027014Gd-1554.4100E − 114.0797E − 1290.751.010271.010270Table 2
Comparison between V.S.O.P. built-in yields and ENDF/B-VII.1 inkeff prediction (CFYs).
IndexFission productFission yieldskeffV.S.O.P.ENDF/B-VII.1Relative difference to ENDF/B-VII.1 (%)V.S.O.P.ENDF/B-VII.1Difference to ENDF/B-VII.1 (pcm)1Xe-1356.6023E − 026.5385E − 020.971.010271.01044−172FP-449.4760E − 019.4760E − 010.001.010271.0102703Xe-1366.2701E − 026.3127E − 02−0.681.010271.0102704Kr-835.3076E − 035.3620E − 03−1.021.010271.0102705Zr-956.4678E − 026.5027E − 02−0.541.010271.0102706Mo-951.6410E − 066.5029E − 02−3.96×10101.010271.009111167Mo-975.9600E − 025.9968E − 02−0.621.010271.0102708Tc-996.1284E − 026.1087E − 020.321.010271.01028−19Ru-1015.0501E − 025.1725E − 02−2.421.010271.01026110Ru-1033.1411E − 023.0309E − 023.511.010271.01042−1511Rh-1051.0199E − 029.6416E − 035.471.010271.01030−312Pd-1087.1032E − 045.4125E − 0423.801.010271.01028−113Ag-1092.9903E − 043.1221E − 04−4.411.010271.01027014Cd-1131.2425E − 041.4038E − 04−12.981.010271.01027015I-1312.8325E − 022.8907E − 02−2.051.010271.01022516Xe-1336.7859E − 026.6991E − 021.281.010271.01032−517Pr-1415.8929E − 025.8470E − 020.781.010271.01028−118Pr-1435.9710E − 025.9558E − 020.251.010271.01029−219Nd-1445.4523E − 025.4996E − 02−0.871.010271.01027020Nd-1453.9339E − 023.9334E − 020.011.010271.01027021Nd-1462.9912E − 022.9969E − 02−0.191.010271.01027022Pm-1472.2701E − 022.2467E − 021.031.010271.01035−823Pm-1491.0888E − 021.0816E − 020.591.010271.01031−424Sm-1505.4130E − 062.9998E − 0794.461.010271.01027025Pm-1514.2044E − 034.1877E − 030.401.010271.01028−126Sm-1522.7057E − 032.6691E − 031.351.010271.01029−227Eu-1531.6264E − 031.5828E − 032.681.010271.01029−228Eu-1553.3025E − 043.2136E − 042.691.010271.01028−129Gd-1561.3517E − 041.4853E − 04−9.881.010271.01027030Gd-1576.4651E − 056.1506E − 054.861.010271.010270IFYs are evaluated by subtracting the total contributions of its precursors from experimental measured CFYs. With the improvement of CFYs measurements, the evaluated IFYs become more precise. It could be seen from the table that IFYs in ENDF/B-VII.1 are lower than the built-in fission yields used in V.S.O.P. Except the large discrepancy in the fission yield of Mo-95, all the impacts from replacing fission yields are lower than 20 pcm. The overall impact is 67 pcm (see Table3) when all the yields are replaced without FP-44. The difference is acceptable in effective multiplication factor calculations when substituting built-in V.S.O.P. fission yields with ENDF/B-VII.1 fission yields. The V.S.O.P. burnup model is further used to conduct fission product yields uncertainty propagation as described in Section 2.3.Table 3
Reference calculation between built-in fission yields of V.S.O.P. and ENDF/B-VII.1 (all substitution without FP-44).
keffDifference to ENDF/B-VII.1 (pcm)Built-in V.S.O.P. fission yieldsENDF/B-VII.1 fission product yields (without update)1.010271.0096067


## 2.3. Uncertainty Quantification Scheme

The HTR-PM [28] reactor core is modelled in V.S.O.P. computer code system to analyse the uncertainty propagation of fission yields in this work. 15 times recirculation of fuel is adopted and the average discharge burnup value is around 90,210MW⋅d/tU with fresh fuel having 8.5 wt.% enrichment. As fuel recirculation tightly couples the neutronics and burnup calculation spatially inside the core, it is difficult to separate the uncertainty propagation step by step. Stochastic UQ method is used to investigate the uncertainty propagation in equilibrium core state. An uncertainty propagation scheme is proposed in this work (see Figure 4).Figure 4
Flow chart of fission product yields uncertainty propagation.Two sampling methods are implemented in this work, namely, normal-based sampling and lognormal-based sampling. Different from normal-based sampling, lognormal-based sampling requires a lognormal transformation of the original mean vector and covariance matrix. When the IFYs samples are generated, their corresponding CFYs are calculated and combine them to form self-consistent fission yield samples. These combined IFYs and CFYs samples are propagated to V.S.O.P. HTR-PM model for further uncertainty analysis. Detailed Bayesian updating method description and lognormal-based sampling procedures will be introduced in Section3.


## 3. Fission Product Yield Perturbation

### 3.1. U-235 Thermal Neutron-Induced IFYs Covariances Estimation

Bayesian updating method or the generalized least square method (GLSM) is a data adjustment method, which allows the prior data being updated with combination of new knowledge about these data. Such knowledge could be measured data or physical constraints imposed on these prior data. The present work applies Bayesian updating method to estimate the covariance matrix of IFYs based on ENDF/B-VII.1 thermal neutron induced U-235 fission yields sublibrary. The specification of this method is briefly recalled as follows.Consider a multivariate linear regression model shown in(2)c=Xy+ε,where c and y∈ℝn×1 are observables and parameters to be updated or estimated, respectively. X∈ℝn×n is the design matrix that represents linear mapping between estimating parameters and observables. ε∈ℝn×1 are the measurement errors of observables with expectation Eε=0∈ℝn×1 and variance Varε=V∈ℝn×n. By the principle of maximum information entropy, it is objective and plausible to assign Gaussian distribution on this error. Similarly, estimating parameters y could also be assigned Gaussian distribution given their expectation Ey=y0 and variance Vary=Z0. The generalized least square problem [29] is formulated by the following minimization in the domain of estimating parameters to find the best least square estimated parameters as(3)minyχ2=Xy−cTV−1Xy−c+y−y0TZ0−1y−y0.The above minimization process could also be interpreted in the perspective of Bayesian updating. Consider the estimated parameters have a prior of Gaussian distribution with densitypy in(4)py=12πn/2Z01/2exp−12y−y0TZ0−1y−y0.And likelihood function determines the probability of any candidate estimated parameters appearing in the observables distribution. Then, likelihood functionpcy is given as(5)pcy=12πn/2V1/2exp−12Xy−cTV−1Xy−c.The posterior distribution of estimated parametersy is therefore calculated by Bayesian theorem and it gives(6)pyc=pcypcpy∝pcypy.Considering the conjugacy between Gaussian prior and likelihood, the posterior estimated parameters follows Gaussian distribution as well. Under quadratic loss, the optimal estimates of true values and their uncertainty are the mean vector and covariance matrix of the posterior distribution. It is worthwhile to mention that the estimated mean vector could maximize the exponential term in (3) and this could also lead to the solution of GLSM in (3).The posterior estimated parameters are obtained as(7)ypost=y+Z1XTV−1c−Xy0,where Z1 is the posterior covariance matrix of estimated parameters and it is shown in (8), and after applying Woodbury matrix identity, it is reformed as (9):(8)Z1=Z0−1+XTV−1X−1,(9)Z1=Z0−Z0XTV+XZ0XT−1XZ0.Here, regarding IFYs as estimated parameters y with prior covariance matrix Z0 (diagonal matrix with only consideration of each fission yields uncertainty in ENDF/B-VII.1), observables c represent the evaluated CFYs in ENDF/B-VII.1, total independent yields, fission system total mass number, and charge number, respectively. The corresponding design matrix could be formulated as follows:(1)
Consistency with CFYs:c=My, where M is the Q-matrix proposed in [11]. It could be formulated from the linear mapping in (1) with the provided branching ratios data in ENDF/B-VII.1 decay sublibrary. This updating process follows Luca Fiorito’s updating procedures [15] on CFYs consistency in JEFF-3.1.2. Different than in previous work [14], this work explicitly constructs this design matrix with branching ratios rather than obtaining each element via direct perturbations using a burnup code. Such procedures allow direct examination of consistency between IFYs and CFYs in the current releases of ENDF/B-VII.1. Total IFYs, total mass number, and total charge number conservations are implemented following the procedures proposed in Pigni et al.’s work [14]. The updating results of IFYs’ covariance matrix are in(10)Z1=Z0−Z0MTV+MZ0MT−1MZ0.(2)
Conservation of binary fission:Ty=UTy, where U∈ℝn×1 is a unity vector. The sum of total yield Ty is 2.0 with summation precision of σsum2=1.0×10−5. The updated covariance matrix subsequent to (10) is listed in (11). It should be noticed that ternary fissions may occur; however they are not considered in ENDF/B-VII.1 and these ternary fissions are not included in this updating process:(11)Z2=Z1−Z1Uσsum2+UTZ1U−1UTZ1.(3)
Conservation of fission system mass number:TM=NTY, where N∈ℝn×1 whose element corresponds to the mass number of each fission product. The total mass number of fission system is conserved to 233.57915 (considering the average prompt neutrons released at 0.0253 eV is 2.42085 recorded in ENDF/B-VII.1 and mass defect of U-235 is not considered). The assumed variance of total mass number is 1.0×10−5. The updated covariance matrix subsequent to (11) is shown in(12)Z3=Z2−Z2Nσsum2+NTZ2N−1NTZ2.(4)
Conservation of fission system charge number:TC=WTy, where W∈ℝn×1 with each element being the charge number of each fission product considered. The total charge number of fission system is conserved as 92.05318. This total charge number is calculated from the charge numbers of each fission product weighted by their corresponding IFYs provided in ENDF/B-VII.1. It is observed in this work that if we take the total charge number as exactly 92.0, the calculated CFYs calculated from updated IFYs will have large discrepancy with CFYs provided in the library. And this discrepancy will be narrowed when we take the decimal digits into consideration. The updated covariance subsequent to (12) is shown in(13)Z4=Z3−Z3Wσsum2+WTZ3W−1WTZ3.Correlation matrix of updated IFYs is plotted (see Figure5). These correlations are introduced sequentially to cooperate the consistency with CFYs, conservation of binary fission, mass number, and charge number of fission system. Figure 5(a) shows that there is a significantly two-humped tendency in the correlation distribution. This tendency is similar with the two-humped distribution of IFYs, where many correlations are introduced from the conservation constraints in fission system, while fewer correlations are introduced between humped part and valley part. And Figure 5(b) presents a close look of the correlations among fission product index range from 65 to 245. It could be noticed that the diagonal of this correlation matrix is divided into several small groups regarding different decay chains. IFYs within each decay chain have negative correlation with each other introduced from the consistency of CFYs.Figure 5
Estimated IFYs correlation matrix. (a) All the estimated correlation information. (b) Section of the estimated correlation information. Red dot indicates the positive correlation and blue dot indicates the negative correlation. The fission product index refers to each fission product identified by its charge numberZ, mass number A, and isomeric state IZZAAAI. These indices are grouped by the mass number and arranged in a descending manner. For each mass group, charge number is ordered in an ascending manner to cooperate the β− decay.
(a)(b)The updated IFYs are compared with the prior fission yields recorded in ENDF/B-VII.1 (see Figure6). It could be seen that small adjustment is introduced to fission product yields in the two-humped part, while larger adjustment is introduced in the valley and two tail parts. This is mainly because IFYs in those parts have smaller prior fission yields and they are not as accurately evaluated as those larger ones in the two-humped part. Therefore, more adjustments are expected in those regions. The updated and prior standard deviations are presented and compared (see Figure 7). It could be seen that the adopted updating procedures could reduce the uncertainty of updated IFYs. This is mainly due to the introduced constraints that further constrain the uncertainty of these fission yields and introduce covariances among them.Figure 6
IFYs distribution before and after sequential Bayesian updating. Natural logarithm is presented on the y-scale. The fission product index refers to each fission product identified by its charge numberZ, mass number A, and isomeric state IZZAAAI. These indices are grouped by the mass number and arranged in a descending manner.Figure 7
IFYs standard deviation (STD) distribution before and after sequential Bayesian updating. Natural logarithm is presented on they-scale. The fission product index refers to each fission product identified by its charge number Z, mass number A, and isomeric state IZZAAAI. These indices are grouped by the mass number and arranged in a descending manner.The final updated covariance matrix of IFYsZ4 and the posterior IFYs mean vector y4 are applied to generate the perturbation samples of IFYs. The detailed sampling procedures are further discussed in the following section.


### 3.2. Lognormal-Based Sampling Procedures

Considering IFYs are inherently positive, random sampling under normal distribution could draw unphysical negative samples. These negative samples would appear significantly when the sampled parameters have large uncertainty (e.g., relative differenceσ/μ>30%). Smith et al. concluded that when the relative uncertainty of a random variable exceeds 30%, the probability distribution of this parameter chosen to represent its physical uncertain information tends to be skewed noticeably [18] and the drawn negative samples fraction tends to grow. It could therefore be concluded that normal distribution is not adequate to describe inherently positive random variables whose uncertainties are large, because it could not capture the skewness of random variable distribution. By the principle of maximum information entropy, lognormal distribution is suggested to be the optimal choice for inherently positive parameter when only expectation and variance are known about this parameter [10, 29]. Larger relative uncertainty would result in a more skewed distribution (shown in Figure 8). Lognormal distribution is shifting to a normal-like distribution as its relative uncertainty becomes lower than 30%, where skewness of the distribution is not significant.Figure 8
Lognormal distribution of random variableX in terms of its relative uncertainty. Relative uncertainty R=σX/μX is ranged from 10% to 80% and μX=2.0. Dashed line shows the distribution with relative uncertainty lower than or equal to 40%, whereas solid line indicates the distribution with relative uncertainty larger than 40%.The updated posterior IFYs relative uncertainties are compared with prior relative uncertainties (see Figure9) in our previous work [30]. Except for a few fission products which have their relative uncertainties increased, most fission products have their corresponding relative uncertainties decreased to around 42%. The increased relative uncertainty fission products are Ag130m0, Cd129m0, Sn127m1, Cd126m0, In126m0, Sb124m1, Zn123m1, Ag115m0, Y93m1, Y93m0, Se85m1, and Ge77m0. Their relative uncertainties increased due to their updated smaller posterior mean values. From Figure 10, it could be seen that most fission yields standard deviations have been reduced because of the updating process. However, the above fission products have their mean value updated even smaller and that makes their relative uncertainties increased. Compared with the listed monitor fission products for fission of U-235 in Fiorito et al.’s work [15], they are not included and we may think they are less relevant to the reactor burnup and critical calculation. When applying simple random sampling procedures under normal distribution, drawing samples in ℝn×S from the Ny4,Z4, where n is the number of fission yields and S is the sample size, it is almost impossible to draw a sample set with all positive yields as the yields domain is too large (e.g., n>900).Figure 9
Relative uncertainties of prior and posterior independent yields in ENDF/B-VII.1 of fission products. These fission products are grouped with their corresponding mass number and the first mass numbers are labelled in this figure [30]. The increased relative uncertainty fission products are Ag130m0, Cd129m0, Sn127m1, Cd126m0, In126m0, Sb124m1, Zn123m1, Ag115m0, Y93m1, Y93m0, Se85m1, and Ge77m0.Figure 10
Posterior to prior fission yields standard deviation ratios. These fission products are grouped with their corresponding mass number and the first mass numbers are labelled in this figure.In this work, lognormal random sampling procedures are applied to generate IFYs perturbation samples. The sampling follows the development of Žerovnik et al. [19] and applies it into the generation of IFYs samples. Multivariate lognormal distribution is defined as(14)L=lny∼Nμl,Zl,where y is the posterior IFYs with expectation y4 and covariance matrix Z4 estimated by Bayesian updating method discussed in Section 3.2 and L∈ℝn×1 is the natural logarithmic value of independent yields. μl and Zl are the corresponding mean and covariance matrix of IFYs in the natural logarithmic domain. The detailed derivation of their relation with parameters in original domain (y4 and Z4) could be found in [20]. The basic idea is recapped in the following.Consider the preservation of probability; the relation between random variables in original domain and logarithmic domain is formulated in(15)pLldl=pYydy.The lognormal distribution density is therefore derived as in(16)pYy=12πn/2Zl1/2∏i=1nyiexp−12lny−μlTZl−1lny−μl.With the logarithmic density function, each element inμl and Zl is derived as(17)μlnyi=lnyi−Varlnyi2.0,(18)covlnxi,lnxj=lncovyi,yjμyiμyj+1,where covyi,yj and μyi are retrieved from the posterior updated IFYs covariance matrix Z4 and updated IFYs mean vector y4. With the calculated distribution parameters μl and Zl, the lognormal-based IFYs sampling procedures could be conducted as follows:(1)
Obtain prior IFYs information including IFYs valuey0 as well as its covariance matrix Z0 from ENDF/B-VII.1 fission yield sublibrary. Implement Bayesian updating procedures detailed in Section 3.2 on the prior information and obtaining the updated IFYs mean vector y4 and the estimated covariance matrix Z4.(2)
Consider IFYs follow lognormal distribution, and transformy4 and Z4 into natural logarithmic domain with (17) and (18). The normal distribution parameters of natural logarithmic yields are obtained as mean vector μl and covariance Zl.(3)
The transformed logarithmic domain covariance could not remain symmetric positive definite (SPD) due to the numerical error in the transformation procedure. A nearest-SPD searching algorithm [31] is therefore applied to search for the nearest SPD approximation of the calculated covariance matrix in the sense of least Frobenius norm difference. The approximated SPD logarithmic domain covariance matrix is thus obtained as Zl′.(4)
Implement the simple random sampling procedures in the logarithmic yield domain with distribution parameters meanμl and approximated SPD covariance matrix Zl′. And the generated logarithmic fission yields sample matrix Pn×S is obtained where n denoted the number of fission products considered and S is sample size.(5)
Take the exponential transformation of each element in sample matrixPn×s and the sampled negative-free samples are generated and denoted as Yn×S.The nearest-SPD searching algorithm approximates non-SPD covariance matrixZl by an approximated matrix Zl′ with relative difference in Frobenius norm Zl−Zl′F/Zl=9.74% and their corresponding eigenvalues distributions are presented in Figure 11. The nearest-SPD searching algorithm could approximate a non-SPD covariance matrix while most of its eigenvalue unchanged.Figure 11
Eigenvalue distribution of relative covariance matrix. Blue dots show the eigenvalue distribution of transformed logarithmic relative covariance. Red dots show the eigenvalue distribution of approximated relative covariance matrix searched by nearest-SPD algorithm. Grey lines show the relative difference of these eigenvalues.The approximation that resides in the above sampling procedures is the SPD approximation of calculated covariance matrix. This approximation could affect consistency of each drawn IFYs sample with the physical constraints imposed on it. There are 1,000 IFYs samples drawn with the lognormal sampling procedures. And the sample mean and standard deviation (STD) for each fission product yield and Pearson’s correlation coefficient between these fission yields are calculated and justified by comparison with its corresponding population values in updatedy4 and Z4 (see Table 4).Table 4
Independent yield samples justification.
Relative differenceMeanRMSEMax of absoluteSample mean8.244×10−41.604×10−21.456×10−1Sample STD6.492×10−36.016×10−28.926×10−1Sampleρ−9.494×10−41.212×10−47.933×100CommentSTD: sample standard deviation;ρ: Pearson’s correlation coefficientRMSE: root mean square errorMax of absolute: the maximum absolute value of relative differenceTable4 indicates that the proposed lognormal sampling procedures could obtain an overall representation of IFYs population distribution considering the lower RMSE. However, there still exist a few fission products listed in Figure 12 having large biases compared with their corresponding population values considering the maximum of absolute relative difference. After comparing these fission products with the monitor fission products for thermal neutron induced U-235 fission listed in Fiorito et al.’s work [15], they are not included and could be considered less relevant to reactor burnup and criticality calculations. These outliers’ appearance could result from the nearest-SPD procedures and a further investigation regarding this will be conducted in future work. Figure 13 presents the sampled Pearson’s correlation coefficients relative difference to their corresponding population values. It could be seen that simple random sampling procedure is not an efficient sampler for sampling low correlation fission yields ρ<0.1 as shown in the neighbour around 0.00 in this figure. However, these low correlations could have little impact on the uncertainty quantification of fission yields compared with large correlations. As for the larger correlations ρ>0.25, 1,000 samples are sufficient for maintaining the Bayesian updated correlations and this discrepancy could be further reduced when increasing the sample size. A more efficient sampler, like Latin Hypercubic Sampler (LHS), could be adopted to guarantee more precise results when using 1,000 samples and this will be adopted in future work.Figure 12
Relative difference of lognormal-based sampled independent fission yields mean values (blue dots) and STD (orange dots) to Bayesian updated values. The outlier fission products are (mean values) Ag130m0, Sn127m1, Cd126m0, In126m0, In118m1, Br86m1, and Ge75m0 and (STD values) Sb131m0, Ag130m0, Cd126m0, Sn122m0, and Br86m1.Figure 13
Relative difference of lognormal-based sampled independent fission yields Pearson’s correlation coefficients to Bayesian updated values against Bayesian updated correlation coefficients. Blue dots represent the values obtained from 1,000 samples and orange dots represent the values obtained from 10,000 samples.The consistency of IFYs samples with these imposed physical constraints is justified in Table5. The conservation parameters (e.g., total fission yields, total mass number, and total charge number) are calculated for each yield samples and the mean and standard deviation are summarized to compare with the target conservation value. It is found that although the consistency is not strictly restored as the standard deviation of the total yield is larger than the imposed 10−5, their mean values are close enough to the target value indicating the constraints are maintained. The large standard deviation is originated from the approximation mentioned above.Table 5
Independent yield physical constraints consistency.
ConstraintTargetMeanSTDBinary fission2.000002.000624.1960×10−3Mass number233.57915233.648664.6685×10−1Charge number92.0531892.076471.8432×10−1In order to examine the difference between normal-based sampling and lognormal-based sampling, 1,000 samples are drawn from the IFYs distribution of Zr95m0, Mo95m0, and Cs134m0. Notationm0 indicates these fission products are at ground state. The IFYs of these three fission products are explicitly involved in V.S.O.P. burnup calculation and are important for reactor decay heat release calculations. Especially for Cs134m0, it is one of the main decay heat contributors of UOX fuels in long-term after reactor shutdown [32]. The updated relative uncertainty of Zr95m0 IFY is 16.1% while Cs134m0 and Mo95m0 have their relative uncertainties of 38.4% and 65.7%, respectively. From the sampled histogram of these fission products IFYs samples (see Figures 14–16), lognormal-based sampling procedures (blue bars) could effectively capture the skewness of these fission yields and permit “negative-free” samples. It is also observed that the skewness of these fission products would become larger as their relative uncertainties become larger (e.g., Mo95m0 and Cs134m0).Figure 14
Histogram comparison of lognormal sampled (a) and normal sampled (b) independent fission yields of Zr95m0. This fission product has relative uncertainty of 16.1%. The text presents the sampled mean value and STD and the population mean and STD for Zr95m0 are1.271856×10−3 and 2.029263×10−4.
(a)(b)Figure 15
Histogram comparison of lognormal sampled (a) and normal sampled (b) independent fission yields of Cs134m0. This fission product has relative uncertainty of 38.4%. The text presents the sampled mean value and STD and the population mean and STD for Cs134m0 are3.8544056×10−8 and 1.456322×10−8.
(a)(b)Figure 16
Histogram comparison of lognormal sampled (a) and normal sampled (b) independent fission yields samples of Mo95m0. This fission product has relative uncertainty of 65.7%. The text presents the sampled mean value and STD and the population mean and STD for Mo95m0 are4.939650×10−12 and 3.161380×10−12.
(a)(b)


## 3.1. U-235 Thermal Neutron-Induced IFYs Covariances Estimation

Bayesian updating method or the generalized least square method (GLSM) is a data adjustment method, which allows the prior data being updated with combination of new knowledge about these data. Such knowledge could be measured data or physical constraints imposed on these prior data. The present work applies Bayesian updating method to estimate the covariance matrix of IFYs based on ENDF/B-VII.1 thermal neutron induced U-235 fission yields sublibrary. The specification of this method is briefly recalled as follows.Consider a multivariate linear regression model shown in(2)c=Xy+ε,where c and y∈ℝn×1 are observables and parameters to be updated or estimated, respectively. X∈ℝn×n is the design matrix that represents linear mapping between estimating parameters and observables. ε∈ℝn×1 are the measurement errors of observables with expectation Eε=0∈ℝn×1 and variance Varε=V∈ℝn×n. By the principle of maximum information entropy, it is objective and plausible to assign Gaussian distribution on this error. Similarly, estimating parameters y could also be assigned Gaussian distribution given their expectation Ey=y0 and variance Vary=Z0. The generalized least square problem [29] is formulated by the following minimization in the domain of estimating parameters to find the best least square estimated parameters as(3)minyχ2=Xy−cTV−1Xy−c+y−y0TZ0−1y−y0.The above minimization process could also be interpreted in the perspective of Bayesian updating. Consider the estimated parameters have a prior of Gaussian distribution with densitypy in(4)py=12πn/2Z01/2exp−12y−y0TZ0−1y−y0.And likelihood function determines the probability of any candidate estimated parameters appearing in the observables distribution. Then, likelihood functionpcy is given as(5)pcy=12πn/2V1/2exp−12Xy−cTV−1Xy−c.The posterior distribution of estimated parametersy is therefore calculated by Bayesian theorem and it gives(6)pyc=pcypcpy∝pcypy.Considering the conjugacy between Gaussian prior and likelihood, the posterior estimated parameters follows Gaussian distribution as well. Under quadratic loss, the optimal estimates of true values and their uncertainty are the mean vector and covariance matrix of the posterior distribution. It is worthwhile to mention that the estimated mean vector could maximize the exponential term in (3) and this could also lead to the solution of GLSM in (3).The posterior estimated parameters are obtained as(7)ypost=y+Z1XTV−1c−Xy0,where Z1 is the posterior covariance matrix of estimated parameters and it is shown in (8), and after applying Woodbury matrix identity, it is reformed as (9):(8)Z1=Z0−1+XTV−1X−1,(9)Z1=Z0−Z0XTV+XZ0XT−1XZ0.Here, regarding IFYs as estimated parameters y with prior covariance matrix Z0 (diagonal matrix with only consideration of each fission yields uncertainty in ENDF/B-VII.1), observables c represent the evaluated CFYs in ENDF/B-VII.1, total independent yields, fission system total mass number, and charge number, respectively. The corresponding design matrix could be formulated as follows:(1)
Consistency with CFYs:c=My, where M is the Q-matrix proposed in [11]. It could be formulated from the linear mapping in (1) with the provided branching ratios data in ENDF/B-VII.1 decay sublibrary. This updating process follows Luca Fiorito’s updating procedures [15] on CFYs consistency in JEFF-3.1.2. Different than in previous work [14], this work explicitly constructs this design matrix with branching ratios rather than obtaining each element via direct perturbations using a burnup code. Such procedures allow direct examination of consistency between IFYs and CFYs in the current releases of ENDF/B-VII.1. Total IFYs, total mass number, and total charge number conservations are implemented following the procedures proposed in Pigni et al.’s work [14]. The updating results of IFYs’ covariance matrix are in(10)Z1=Z0−Z0MTV+MZ0MT−1MZ0.(2)
Conservation of binary fission:Ty=UTy, where U∈ℝn×1 is a unity vector. The sum of total yield Ty is 2.0 with summation precision of σsum2=1.0×10−5. The updated covariance matrix subsequent to (10) is listed in (11). It should be noticed that ternary fissions may occur; however they are not considered in ENDF/B-VII.1 and these ternary fissions are not included in this updating process:(11)Z2=Z1−Z1Uσsum2+UTZ1U−1UTZ1.(3)
Conservation of fission system mass number:TM=NTY, where N∈ℝn×1 whose element corresponds to the mass number of each fission product. The total mass number of fission system is conserved to 233.57915 (considering the average prompt neutrons released at 0.0253 eV is 2.42085 recorded in ENDF/B-VII.1 and mass defect of U-235 is not considered). The assumed variance of total mass number is 1.0×10−5. The updated covariance matrix subsequent to (11) is shown in(12)Z3=Z2−Z2Nσsum2+NTZ2N−1NTZ2.(4)
Conservation of fission system charge number:TC=WTy, where W∈ℝn×1 with each element being the charge number of each fission product considered. The total charge number of fission system is conserved as 92.05318. This total charge number is calculated from the charge numbers of each fission product weighted by their corresponding IFYs provided in ENDF/B-VII.1. It is observed in this work that if we take the total charge number as exactly 92.0, the calculated CFYs calculated from updated IFYs will have large discrepancy with CFYs provided in the library. And this discrepancy will be narrowed when we take the decimal digits into consideration. The updated covariance subsequent to (12) is shown in(13)Z4=Z3−Z3Wσsum2+WTZ3W−1WTZ3.Correlation matrix of updated IFYs is plotted (see Figure5). These correlations are introduced sequentially to cooperate the consistency with CFYs, conservation of binary fission, mass number, and charge number of fission system. Figure 5(a) shows that there is a significantly two-humped tendency in the correlation distribution. This tendency is similar with the two-humped distribution of IFYs, where many correlations are introduced from the conservation constraints in fission system, while fewer correlations are introduced between humped part and valley part. And Figure 5(b) presents a close look of the correlations among fission product index range from 65 to 245. It could be noticed that the diagonal of this correlation matrix is divided into several small groups regarding different decay chains. IFYs within each decay chain have negative correlation with each other introduced from the consistency of CFYs.Figure 5
Estimated IFYs correlation matrix. (a) All the estimated correlation information. (b) Section of the estimated correlation information. Red dot indicates the positive correlation and blue dot indicates the negative correlation. The fission product index refers to each fission product identified by its charge numberZ, mass number A, and isomeric state IZZAAAI. These indices are grouped by the mass number and arranged in a descending manner. For each mass group, charge number is ordered in an ascending manner to cooperate the β− decay.
(a)(b)The updated IFYs are compared with the prior fission yields recorded in ENDF/B-VII.1 (see Figure6). It could be seen that small adjustment is introduced to fission product yields in the two-humped part, while larger adjustment is introduced in the valley and two tail parts. This is mainly because IFYs in those parts have smaller prior fission yields and they are not as accurately evaluated as those larger ones in the two-humped part. Therefore, more adjustments are expected in those regions. The updated and prior standard deviations are presented and compared (see Figure 7). It could be seen that the adopted updating procedures could reduce the uncertainty of updated IFYs. This is mainly due to the introduced constraints that further constrain the uncertainty of these fission yields and introduce covariances among them.Figure 6
IFYs distribution before and after sequential Bayesian updating. Natural logarithm is presented on the y-scale. The fission product index refers to each fission product identified by its charge numberZ, mass number A, and isomeric state IZZAAAI. These indices are grouped by the mass number and arranged in a descending manner.Figure 7
IFYs standard deviation (STD) distribution before and after sequential Bayesian updating. Natural logarithm is presented on they-scale. The fission product index refers to each fission product identified by its charge number Z, mass number A, and isomeric state IZZAAAI. These indices are grouped by the mass number and arranged in a descending manner.The final updated covariance matrix of IFYsZ4 and the posterior IFYs mean vector y4 are applied to generate the perturbation samples of IFYs. The detailed sampling procedures are further discussed in the following section.


## 3.2. Lognormal-Based Sampling Procedures

Considering IFYs are inherently positive, random sampling under normal distribution could draw unphysical negative samples. These negative samples would appear significantly when the sampled parameters have large uncertainty (e.g., relative differenceσ/μ>30%). Smith et al. concluded that when the relative uncertainty of a random variable exceeds 30%, the probability distribution of this parameter chosen to represent its physical uncertain information tends to be skewed noticeably [18] and the drawn negative samples fraction tends to grow. It could therefore be concluded that normal distribution is not adequate to describe inherently positive random variables whose uncertainties are large, because it could not capture the skewness of random variable distribution. By the principle of maximum information entropy, lognormal distribution is suggested to be the optimal choice for inherently positive parameter when only expectation and variance are known about this parameter [10, 29]. Larger relative uncertainty would result in a more skewed distribution (shown in Figure 8). Lognormal distribution is shifting to a normal-like distribution as its relative uncertainty becomes lower than 30%, where skewness of the distribution is not significant.Figure 8
Lognormal distribution of random variableX in terms of its relative uncertainty. Relative uncertainty R=σX/μX is ranged from 10% to 80% and μX=2.0. Dashed line shows the distribution with relative uncertainty lower than or equal to 40%, whereas solid line indicates the distribution with relative uncertainty larger than 40%.The updated posterior IFYs relative uncertainties are compared with prior relative uncertainties (see Figure9) in our previous work [30]. Except for a few fission products which have their relative uncertainties increased, most fission products have their corresponding relative uncertainties decreased to around 42%. The increased relative uncertainty fission products are Ag130m0, Cd129m0, Sn127m1, Cd126m0, In126m0, Sb124m1, Zn123m1, Ag115m0, Y93m1, Y93m0, Se85m1, and Ge77m0. Their relative uncertainties increased due to their updated smaller posterior mean values. From Figure 10, it could be seen that most fission yields standard deviations have been reduced because of the updating process. However, the above fission products have their mean value updated even smaller and that makes their relative uncertainties increased. Compared with the listed monitor fission products for fission of U-235 in Fiorito et al.’s work [15], they are not included and we may think they are less relevant to the reactor burnup and critical calculation. When applying simple random sampling procedures under normal distribution, drawing samples in ℝn×S from the Ny4,Z4, where n is the number of fission yields and S is the sample size, it is almost impossible to draw a sample set with all positive yields as the yields domain is too large (e.g., n>900).Figure 9
Relative uncertainties of prior and posterior independent yields in ENDF/B-VII.1 of fission products. These fission products are grouped with their corresponding mass number and the first mass numbers are labelled in this figure [30]. The increased relative uncertainty fission products are Ag130m0, Cd129m0, Sn127m1, Cd126m0, In126m0, Sb124m1, Zn123m1, Ag115m0, Y93m1, Y93m0, Se85m1, and Ge77m0.Figure 10
Posterior to prior fission yields standard deviation ratios. These fission products are grouped with their corresponding mass number and the first mass numbers are labelled in this figure.In this work, lognormal random sampling procedures are applied to generate IFYs perturbation samples. The sampling follows the development of Žerovnik et al. [19] and applies it into the generation of IFYs samples. Multivariate lognormal distribution is defined as(14)L=lny∼Nμl,Zl,where y is the posterior IFYs with expectation y4 and covariance matrix Z4 estimated by Bayesian updating method discussed in Section 3.2 and L∈ℝn×1 is the natural logarithmic value of independent yields. μl and Zl are the corresponding mean and covariance matrix of IFYs in the natural logarithmic domain. The detailed derivation of their relation with parameters in original domain (y4 and Z4) could be found in [20]. The basic idea is recapped in the following.Consider the preservation of probability; the relation between random variables in original domain and logarithmic domain is formulated in(15)pLldl=pYydy.The lognormal distribution density is therefore derived as in(16)pYy=12πn/2Zl1/2∏i=1nyiexp−12lny−μlTZl−1lny−μl.With the logarithmic density function, each element inμl and Zl is derived as(17)μlnyi=lnyi−Varlnyi2.0,(18)covlnxi,lnxj=lncovyi,yjμyiμyj+1,where covyi,yj and μyi are retrieved from the posterior updated IFYs covariance matrix Z4 and updated IFYs mean vector y4. With the calculated distribution parameters μl and Zl, the lognormal-based IFYs sampling procedures could be conducted as follows:(1)
Obtain prior IFYs information including IFYs valuey0 as well as its covariance matrix Z0 from ENDF/B-VII.1 fission yield sublibrary. Implement Bayesian updating procedures detailed in Section 3.2 on the prior information and obtaining the updated IFYs mean vector y4 and the estimated covariance matrix Z4.(2)
Consider IFYs follow lognormal distribution, and transformy4 and Z4 into natural logarithmic domain with (17) and (18). The normal distribution parameters of natural logarithmic yields are obtained as mean vector μl and covariance Zl.(3)
The transformed logarithmic domain covariance could not remain symmetric positive definite (SPD) due to the numerical error in the transformation procedure. A nearest-SPD searching algorithm [31] is therefore applied to search for the nearest SPD approximation of the calculated covariance matrix in the sense of least Frobenius norm difference. The approximated SPD logarithmic domain covariance matrix is thus obtained as Zl′.(4)
Implement the simple random sampling procedures in the logarithmic yield domain with distribution parameters meanμl and approximated SPD covariance matrix Zl′. And the generated logarithmic fission yields sample matrix Pn×S is obtained where n denoted the number of fission products considered and S is sample size.(5)
Take the exponential transformation of each element in sample matrixPn×s and the sampled negative-free samples are generated and denoted as Yn×S.The nearest-SPD searching algorithm approximates non-SPD covariance matrixZl by an approximated matrix Zl′ with relative difference in Frobenius norm Zl−Zl′F/Zl=9.74% and their corresponding eigenvalues distributions are presented in Figure 11. The nearest-SPD searching algorithm could approximate a non-SPD covariance matrix while most of its eigenvalue unchanged.Figure 11
Eigenvalue distribution of relative covariance matrix. Blue dots show the eigenvalue distribution of transformed logarithmic relative covariance. Red dots show the eigenvalue distribution of approximated relative covariance matrix searched by nearest-SPD algorithm. Grey lines show the relative difference of these eigenvalues.The approximation that resides in the above sampling procedures is the SPD approximation of calculated covariance matrix. This approximation could affect consistency of each drawn IFYs sample with the physical constraints imposed on it. There are 1,000 IFYs samples drawn with the lognormal sampling procedures. And the sample mean and standard deviation (STD) for each fission product yield and Pearson’s correlation coefficient between these fission yields are calculated and justified by comparison with its corresponding population values in updatedy4 and Z4 (see Table 4).Table 4
Independent yield samples justification.
Relative differenceMeanRMSEMax of absoluteSample mean8.244×10−41.604×10−21.456×10−1Sample STD6.492×10−36.016×10−28.926×10−1Sampleρ−9.494×10−41.212×10−47.933×100CommentSTD: sample standard deviation;ρ: Pearson’s correlation coefficientRMSE: root mean square errorMax of absolute: the maximum absolute value of relative differenceTable4 indicates that the proposed lognormal sampling procedures could obtain an overall representation of IFYs population distribution considering the lower RMSE. However, there still exist a few fission products listed in Figure 12 having large biases compared with their corresponding population values considering the maximum of absolute relative difference. After comparing these fission products with the monitor fission products for thermal neutron induced U-235 fission listed in Fiorito et al.’s work [15], they are not included and could be considered less relevant to reactor burnup and criticality calculations. These outliers’ appearance could result from the nearest-SPD procedures and a further investigation regarding this will be conducted in future work. Figure 13 presents the sampled Pearson’s correlation coefficients relative difference to their corresponding population values. It could be seen that simple random sampling procedure is not an efficient sampler for sampling low correlation fission yields ρ<0.1 as shown in the neighbour around 0.00 in this figure. However, these low correlations could have little impact on the uncertainty quantification of fission yields compared with large correlations. As for the larger correlations ρ>0.25, 1,000 samples are sufficient for maintaining the Bayesian updated correlations and this discrepancy could be further reduced when increasing the sample size. A more efficient sampler, like Latin Hypercubic Sampler (LHS), could be adopted to guarantee more precise results when using 1,000 samples and this will be adopted in future work.Figure 12
Relative difference of lognormal-based sampled independent fission yields mean values (blue dots) and STD (orange dots) to Bayesian updated values. The outlier fission products are (mean values) Ag130m0, Sn127m1, Cd126m0, In126m0, In118m1, Br86m1, and Ge75m0 and (STD values) Sb131m0, Ag130m0, Cd126m0, Sn122m0, and Br86m1.Figure 13
Relative difference of lognormal-based sampled independent fission yields Pearson’s correlation coefficients to Bayesian updated values against Bayesian updated correlation coefficients. Blue dots represent the values obtained from 1,000 samples and orange dots represent the values obtained from 10,000 samples.The consistency of IFYs samples with these imposed physical constraints is justified in Table5. The conservation parameters (e.g., total fission yields, total mass number, and total charge number) are calculated for each yield samples and the mean and standard deviation are summarized to compare with the target conservation value. It is found that although the consistency is not strictly restored as the standard deviation of the total yield is larger than the imposed 10−5, their mean values are close enough to the target value indicating the constraints are maintained. The large standard deviation is originated from the approximation mentioned above.Table 5
Independent yield physical constraints consistency.
ConstraintTargetMeanSTDBinary fission2.000002.000624.1960×10−3Mass number233.57915233.648664.6685×10−1Charge number92.0531892.076471.8432×10−1In order to examine the difference between normal-based sampling and lognormal-based sampling, 1,000 samples are drawn from the IFYs distribution of Zr95m0, Mo95m0, and Cs134m0. Notationm0 indicates these fission products are at ground state. The IFYs of these three fission products are explicitly involved in V.S.O.P. burnup calculation and are important for reactor decay heat release calculations. Especially for Cs134m0, it is one of the main decay heat contributors of UOX fuels in long-term after reactor shutdown [32]. The updated relative uncertainty of Zr95m0 IFY is 16.1% while Cs134m0 and Mo95m0 have their relative uncertainties of 38.4% and 65.7%, respectively. From the sampled histogram of these fission products IFYs samples (see Figures 14–16), lognormal-based sampling procedures (blue bars) could effectively capture the skewness of these fission yields and permit “negative-free” samples. It is also observed that the skewness of these fission products would become larger as their relative uncertainties become larger (e.g., Mo95m0 and Cs134m0).Figure 14
Histogram comparison of lognormal sampled (a) and normal sampled (b) independent fission yields of Zr95m0. This fission product has relative uncertainty of 16.1%. The text presents the sampled mean value and STD and the population mean and STD for Zr95m0 are1.271856×10−3 and 2.029263×10−4.
(a)(b)Figure 15
Histogram comparison of lognormal sampled (a) and normal sampled (b) independent fission yields of Cs134m0. This fission product has relative uncertainty of 38.4%. The text presents the sampled mean value and STD and the population mean and STD for Cs134m0 are3.8544056×10−8 and 1.456322×10−8.
(a)(b)Figure 16
Histogram comparison of lognormal sampled (a) and normal sampled (b) independent fission yields samples of Mo95m0. This fission product has relative uncertainty of 65.7%. The text presents the sampled mean value and STD and the population mean and STD for Mo95m0 are4.939650×10−12 and 3.161380×10−12.
(a)(b)


## 4. Results and Discussion

### 4.1. Uncertainty Analysis of the Effective Multiplication Factor at Equilibrium Core

The unperturbed burnup calculation is conducted with V.S.O.P. built-in fission yields library and ENDF/B-VII.1 posterior fission yields. Figure17 shows that reactor achieved the equilibrium state after operating longer than 2500 days. Effective multiplication factor calculated from ENDF/B-VII.1 posterior fission yields is compared with that calculated from V.S.O.P. built-in fission yields and the total discrepancy at equilibrium core state (which is at the end point of fuel cycle time in Figure 17) is lower than 50 pcm which is small enough for the following fission product yields uncertainty propagation analysis. The comparison between ENDF/B-VII.1 posterior fission yields predicated keff (black dashed line) and built-in yields predicted keff (orange dashed line) are shown in Figure 18(b). This discrepancy is within the sampling distribution of keff.Figure 17
Multiplication factor predictions with operation time. Operation time involves a running-phase (0–2500 days) and equilibrium core state (3049 days). The multiplication factor prediction differences of V.S.O.P. built-in yields and ENDF/B-VII.1 posterior yields are plotted.Figure 18
Histogram of 1,000 multiplication factor samples obtained at equilibrium core state. (a) Samples obtained by normal-based sampling procedures. (b) Samples obtained by lognormal-based samples procedures. Red line indicates the superimposed fitted normal density from the sampled data.
(a)(b)1,000 fission yields samples are generated with normal-based sampling procedures and lognormal-based sampling procedures and they are propagated to V.S.O.P. burnup calculation to obtainkeff samples under equilibrium core state (3049 days). The sample distributions from these two sampling procedures are drawn and compared (see Figure 18). It is obvious from the comparison that normal-based samples contain fewer distribution information compared with lognormal samples as its distribution range is smaller than that in lognormal samples. This is due to the zero cut-off procedure of the uncontrolled negative samples. Such procedure artificially omits certain information in the original fission yields distributions and could not provide a correspondingly reasonable and satisfied sampling distribution of keff. In this sense, lognormal sampling procedures overcome this problem by imposing a more plausible distribution on fission yields and allow the generation of smaller perturbed samples. Therefore, it leads to a negative skewness (long tail in left) of effective multiplication factor distribution and permits a more rational and persuasive sampling distribution.The uncertainty analysis results are presented (see Table6). The propagated sampled distribution of keff from normal-based sampling method passes the normality test with p value 0.3737 and the quantified relative uncertainty is around 1.09×10−4. Lognormal samples provide a skewed keff distribution and fails the normality test with p value smaller than 0.05. The quantified relative uncertainty from this distribution is 2.58×10−4. The keff quantified from lognormal-based sampling method is larger than that from normal-based sampling method, and this shows that the zero cut-off effect in normal-based sampling method could cause underestimation of fission product yields uncertainty contribution to QoIs.Table 6
Uncertainty analysis results and comparison of multiplication factor at equilibrium core state.
Sampling proceduresNominal prediction V.S.O.P with ENDF/B-VII.1 posterior fission yieldsFission yield uncertainty analysis resultsSamples meanRelative uncertainty95% CINormality test(2)Normal1.001061.001051.09E − 04[1.05E − 04, 1.14E − 04]Passed:p=3.737E−01Lognormal1.001061.001052.44E − 04[2.44E − 04, 2.76E − 04](1)Failed:p=1.103E−04Comment(1) 95% confidence interval is estimated by bootstrapping method with 100,000 bootstrap samples(2) Normality test is conducted on the z-scores of multiplication factor samples with the K-S test


### 4.2. Uncertainty Analysis of Certain Fission Products Atomic Densities

In this section, fission products Zr95m0, Mo95mo, and Cs134m0 atomic densities uncertainties contributed from fission products yields are quantified. Specifically, their uncertainties differences from different sampling methods are compared and discussed. From the discussions in Section3.2, Zr95m0, Cs134m0, and Mo95m0 have their IFYs relative uncertainties of 16.1%, 38.4%, and 65.7%, respectively. And their IFYs sampling results shown in Figures 14–16 (blue bars) indicate that the skewness of the sampled distribution becomes significant with their fission yields relative uncertainties increasing. When we adopt normal-based sampling procedures to a skewed distribution, there will be more negative samples values and the zero cut-off would deliver more underestimation into the uncertainty analysis results.Figures19–21 track the atomic density of Zr95m0, Mo95m0, and Cs134m0 in the loaded 9.8 kg fresh fuel (with enrichment 8.5 wt.%) in HTR-PM along with their irradiation. The horizontal axis indicates the average burnup values of these fuels. As HTR-PM allows recirculation of fuels, 15 times recirculation is adopted in this analysis, which indicates these fresh fuels will be reloaded into the core 15 times before they are finally discharged. The discharged burnup value or the end point of the horizontal axis is 90210.44 MW·d/tU. Throughout the burnup process, the thermal power of reactor core is kept at 250 MW.Figure 19
Relative uncertainties of Zr95m0 against fuels average burnup values. Lognormal-based sampling results (blue line) and normal-based sampling results (orange line) are plotted. The shades in this figure are the 95% confidence interval of relative uncertainties.Figure 20
Relative uncertainties of Mo95m0 against fuels average burnup values. Lognormal-based sampling results (blue line) and normal-based sampling results (orange line) are plotted. The shades in this figure are the 95% confidence interval of relative uncertainties.Figure 21
Relative uncertainties of Cs134m0 against fuels average burnup values. Lognormal-based sampling results (blue line) and normal-based sampling results (orange line) are plotted. The shades in this figure are the 95% confidence interval of relative uncertainties.The atomic densities of Zr-95m0 fluctuate along with the increases of fuels burnup value. This fluctuation is due to the fuel recirculation procedures adopted in V.S.O.P. burnup calculations. There are total 14 lower valleys that appeared in dashed line of Figure19 which corresponds to the 14 times reloading of the fuels from the bottom of the core to the top. For each reloading, the fuels will be irradiated again during their passes through the core. As it could be seen from Figures 19–21, except for the atomic densities’ decrease of Zr95m0 along with the increase of average burnup values, Mo95m0 and Cs134m0 have their atomic densities accumulated throughout the whole burnup process. During the burnup process, their atomic densities relative uncertainties contributed from thermal neutron induced U-235 fission products yields are investigated.The atomic density relative uncertainties of the above three fission products varying with the average burnup values of fuels are plotted in Figures19–21 (blue and orange solid line). The blue and orange shadings in these figures are the 95% confidence interval of relative uncertainty computed by bootstrap method. From these figures, it is worth to mention that lognormal-based sampling quantified atomic density relative uncertainties are larger than that quantified from normal-based sampling for all of these three fission products. This is reasonable as zero cut-off adopted in normal-based sampling method would artificially omit some information provided by fission yields distributions, and this would result in an underestimated atomic density relative uncertainty quantification result. After closely comparing the atomic density relative uncertainties underestimation for Zr95m0 and Cs34m0, it could be seen that this underestimation effect will be enlarged when the fission products IFYs have larger relative uncertainties (Zr95m0: 16.1% and Cs134m0: 38.4%). This is because lognormal distribution would resemble normal distribution when the random variate has smaller relative uncertainty, as discussed in Section 3.2. And in this case, lognormal-based sampling results would be in agreement with those calculated from normal-based sampling. Therefore, this underestimation would be narrowed.Besides, another interesting phenomenon is observed here. This underestimation seems not positively correlated with the relative uncertainty of random variates, as it is seen from comparison between Mo95m0 and Cs134m0. Although Mo95m0 has its relative uncertainty (65.7%) larger than Cs134m0 (38.4%), the underestimation effect observed from Figures20 and 21 shows that the underestimation effect of Mo95m0 is smaller than that of Cs134m0. One possible reason could be the decay of these fission products. As Mo95m0 is the direct descendant of Zr95m0 whose half-life is around 64 days, its atomic density relative uncertainty is contributed both from its own fission yields uncertainty and the atomic density uncertainty of Zr95m0. As Zr95m0 atomic density uncertainty is less underestimated, the atomic density relative uncertainty underestimation in Mo95m0 is therefore counterbalanced. While Cs134m0 is treated as stable fission products in V.S.O.P. burnup fission product chains, its atomic density relative uncertainty is directly related to its fission yields uncertainty and large atomic density relative uncertainty underestimation could be seen. The atomic density relative uncertainties of all these three fission products quantified at 90210.44 MW·d/tU are summarized in Table 7.Table 7
Uncertainty analysis results of Zr95m0, Mo95m0, and Cs134m0.
Fission productsAtomic density (atoms/barn·cm)Fission yield uncertainty analysis resultsNormal-based relative uncertainty (%)Lognormal-based relative uncertainty (%)Relative difference to lognormal-based relative uncertainty (%)Zr95m01.3532×10−70.2420.244−0.8Mo95m04.1478×10−60.4560.458−0.4Cs134m02.3575×10−70.2440.313−22.0CommentThese results show atomic densities in 9.8 kg heavy metal irradiated up to 90210.44 MWd/tU


## 4.1. Uncertainty Analysis of the Effective Multiplication Factor at Equilibrium Core

The unperturbed burnup calculation is conducted with V.S.O.P. built-in fission yields library and ENDF/B-VII.1 posterior fission yields. Figure17 shows that reactor achieved the equilibrium state after operating longer than 2500 days. Effective multiplication factor calculated from ENDF/B-VII.1 posterior fission yields is compared with that calculated from V.S.O.P. built-in fission yields and the total discrepancy at equilibrium core state (which is at the end point of fuel cycle time in Figure 17) is lower than 50 pcm which is small enough for the following fission product yields uncertainty propagation analysis. The comparison between ENDF/B-VII.1 posterior fission yields predicated keff (black dashed line) and built-in yields predicted keff (orange dashed line) are shown in Figure 18(b). This discrepancy is within the sampling distribution of keff.Figure 17
Multiplication factor predictions with operation time. Operation time involves a running-phase (0–2500 days) and equilibrium core state (3049 days). The multiplication factor prediction differences of V.S.O.P. built-in yields and ENDF/B-VII.1 posterior yields are plotted.Figure 18
Histogram of 1,000 multiplication factor samples obtained at equilibrium core state. (a) Samples obtained by normal-based sampling procedures. (b) Samples obtained by lognormal-based samples procedures. Red line indicates the superimposed fitted normal density from the sampled data.
(a)(b)1,000 fission yields samples are generated with normal-based sampling procedures and lognormal-based sampling procedures and they are propagated to V.S.O.P. burnup calculation to obtainkeff samples under equilibrium core state (3049 days). The sample distributions from these two sampling procedures are drawn and compared (see Figure 18). It is obvious from the comparison that normal-based samples contain fewer distribution information compared with lognormal samples as its distribution range is smaller than that in lognormal samples. This is due to the zero cut-off procedure of the uncontrolled negative samples. Such procedure artificially omits certain information in the original fission yields distributions and could not provide a correspondingly reasonable and satisfied sampling distribution of keff. In this sense, lognormal sampling procedures overcome this problem by imposing a more plausible distribution on fission yields and allow the generation of smaller perturbed samples. Therefore, it leads to a negative skewness (long tail in left) of effective multiplication factor distribution and permits a more rational and persuasive sampling distribution.The uncertainty analysis results are presented (see Table6). The propagated sampled distribution of keff from normal-based sampling method passes the normality test with p value 0.3737 and the quantified relative uncertainty is around 1.09×10−4. Lognormal samples provide a skewed keff distribution and fails the normality test with p value smaller than 0.05. The quantified relative uncertainty from this distribution is 2.58×10−4. The keff quantified from lognormal-based sampling method is larger than that from normal-based sampling method, and this shows that the zero cut-off effect in normal-based sampling method could cause underestimation of fission product yields uncertainty contribution to QoIs.Table 6
Uncertainty analysis results and comparison of multiplication factor at equilibrium core state.
Sampling proceduresNominal prediction V.S.O.P with ENDF/B-VII.1 posterior fission yieldsFission yield uncertainty analysis resultsSamples meanRelative uncertainty95% CINormality test(2)Normal1.001061.001051.09E − 04[1.05E − 04, 1.14E − 04]Passed:p=3.737E−01Lognormal1.001061.001052.44E − 04[2.44E − 04, 2.76E − 04](1)Failed:p=1.103E−04Comment(1) 95% confidence interval is estimated by bootstrapping method with 100,000 bootstrap samples(2) Normality test is conducted on the z-scores of multiplication factor samples with the K-S test


## 4.2. Uncertainty Analysis of Certain Fission Products Atomic Densities

In this section, fission products Zr95m0, Mo95mo, and Cs134m0 atomic densities uncertainties contributed from fission products yields are quantified. Specifically, their uncertainties differences from different sampling methods are compared and discussed. From the discussions in Section3.2, Zr95m0, Cs134m0, and Mo95m0 have their IFYs relative uncertainties of 16.1%, 38.4%, and 65.7%, respectively. And their IFYs sampling results shown in Figures 14–16 (blue bars) indicate that the skewness of the sampled distribution becomes significant with their fission yields relative uncertainties increasing. When we adopt normal-based sampling procedures to a skewed distribution, there will be more negative samples values and the zero cut-off would deliver more underestimation into the uncertainty analysis results.Figures19–21 track the atomic density of Zr95m0, Mo95m0, and Cs134m0 in the loaded 9.8 kg fresh fuel (with enrichment 8.5 wt.%) in HTR-PM along with their irradiation. The horizontal axis indicates the average burnup values of these fuels. As HTR-PM allows recirculation of fuels, 15 times recirculation is adopted in this analysis, which indicates these fresh fuels will be reloaded into the core 15 times before they are finally discharged. The discharged burnup value or the end point of the horizontal axis is 90210.44 MW·d/tU. Throughout the burnup process, the thermal power of reactor core is kept at 250 MW.Figure 19
Relative uncertainties of Zr95m0 against fuels average burnup values. Lognormal-based sampling results (blue line) and normal-based sampling results (orange line) are plotted. The shades in this figure are the 95% confidence interval of relative uncertainties.Figure 20
Relative uncertainties of Mo95m0 against fuels average burnup values. Lognormal-based sampling results (blue line) and normal-based sampling results (orange line) are plotted. The shades in this figure are the 95% confidence interval of relative uncertainties.Figure 21
Relative uncertainties of Cs134m0 against fuels average burnup values. Lognormal-based sampling results (blue line) and normal-based sampling results (orange line) are plotted. The shades in this figure are the 95% confidence interval of relative uncertainties.The atomic densities of Zr-95m0 fluctuate along with the increases of fuels burnup value. This fluctuation is due to the fuel recirculation procedures adopted in V.S.O.P. burnup calculations. There are total 14 lower valleys that appeared in dashed line of Figure19 which corresponds to the 14 times reloading of the fuels from the bottom of the core to the top. For each reloading, the fuels will be irradiated again during their passes through the core. As it could be seen from Figures 19–21, except for the atomic densities’ decrease of Zr95m0 along with the increase of average burnup values, Mo95m0 and Cs134m0 have their atomic densities accumulated throughout the whole burnup process. During the burnup process, their atomic densities relative uncertainties contributed from thermal neutron induced U-235 fission products yields are investigated.The atomic density relative uncertainties of the above three fission products varying with the average burnup values of fuels are plotted in Figures19–21 (blue and orange solid line). The blue and orange shadings in these figures are the 95% confidence interval of relative uncertainty computed by bootstrap method. From these figures, it is worth to mention that lognormal-based sampling quantified atomic density relative uncertainties are larger than that quantified from normal-based sampling for all of these three fission products. This is reasonable as zero cut-off adopted in normal-based sampling method would artificially omit some information provided by fission yields distributions, and this would result in an underestimated atomic density relative uncertainty quantification result. After closely comparing the atomic density relative uncertainties underestimation for Zr95m0 and Cs34m0, it could be seen that this underestimation effect will be enlarged when the fission products IFYs have larger relative uncertainties (Zr95m0: 16.1% and Cs134m0: 38.4%). This is because lognormal distribution would resemble normal distribution when the random variate has smaller relative uncertainty, as discussed in Section 3.2. And in this case, lognormal-based sampling results would be in agreement with those calculated from normal-based sampling. Therefore, this underestimation would be narrowed.Besides, another interesting phenomenon is observed here. This underestimation seems not positively correlated with the relative uncertainty of random variates, as it is seen from comparison between Mo95m0 and Cs134m0. Although Mo95m0 has its relative uncertainty (65.7%) larger than Cs134m0 (38.4%), the underestimation effect observed from Figures20 and 21 shows that the underestimation effect of Mo95m0 is smaller than that of Cs134m0. One possible reason could be the decay of these fission products. As Mo95m0 is the direct descendant of Zr95m0 whose half-life is around 64 days, its atomic density relative uncertainty is contributed both from its own fission yields uncertainty and the atomic density uncertainty of Zr95m0. As Zr95m0 atomic density uncertainty is less underestimated, the atomic density relative uncertainty underestimation in Mo95m0 is therefore counterbalanced. While Cs134m0 is treated as stable fission products in V.S.O.P. burnup fission product chains, its atomic density relative uncertainty is directly related to its fission yields uncertainty and large atomic density relative uncertainty underestimation could be seen. The atomic density relative uncertainties of all these three fission products quantified at 90210.44 MW·d/tU are summarized in Table 7.Table 7
Uncertainty analysis results of Zr95m0, Mo95m0, and Cs134m0.
Fission productsAtomic density (atoms/barn·cm)Fission yield uncertainty analysis resultsNormal-based relative uncertainty (%)Lognormal-based relative uncertainty (%)Relative difference to lognormal-based relative uncertainty (%)Zr95m01.3532×10−70.2420.244−0.8Mo95m04.1478×10−60.4560.458−0.4Cs134m02.3575×10−70.2440.313−22.0CommentThese results show atomic densities in 9.8 kg heavy metal irradiated up to 90210.44 MWd/tU


## 5. Conclusions

The present work proposed a stochastic UQ method for propagation fission products yields uncertainties. V.S.O.P. code [26] is used to conduct the burnup calculation of HTR-PM reactor core with allowing 15 times recirculation of fuel pebbles [30]. Uncertainties of thermal neutron induced U-235 IFYs are investigated in this work based on ENDF/B-VII.1. Bayesian updating method is applied to estimate the covariance matrix of IFYs. Lognormal-based sampling method is implemented to generate perturbations of yields samples. The differences of quantified uncertainties between conventional normal-based sampling method and lognormal-based method are addressed and investigated. Specifically, the effect of zero cut-off procedures used in normal-based sampling method is studied and discussed. From the above investigation, conclusions are summarized as follows:(1)
Lognormal-based sampling method could effectively overcome the negative samples generation caused by the large relative uncertainties in fission yields data. Compared with normal-based sampling method, it could provide reasonable and negative-free fission yields samples to permit a more plausible and reasonable QoI sampling distribution for further uncertainty analysis.(2)
The contribution of thermal neutron induced U-235 fission yields uncertainties in ENDF/B-VII.1 tokeff of pebble-bed HTGR at equilibrium core is 0.0258%. This contribution is smaller than that from neutron cross section 0.48% at equilibrium core [33].(3)
The zero cut-off procedures used in conventional normal-based sampling method to overcome the negative fission yields samples appearance would underestimate the uncertainty analysis results. For relative uncertainty of effective multiplication factor, it would underestimate the results by 0.0149% which is around 42% of results obtained from lognormal-based sampling method. For atomic density relative uncertainty, the underestimations are also observed, and especially for Cs134m0, this zero cut-off effect would underestimate the atomic density relative uncertainty by 22% compared with lognormal-based quantified results.It is worth to mention that there are several approximations and simplifications made during the Bayesian updating process and implementing of lognormal-based sampling methods in this work. The considered constraints for Bayesian updating independent yields covariance matrix are preliminary in this work and a more complete and comprehensive study regarding this will be conducted in future work. Also, the effect of using nearest SPD algorithm in implementing lognormal-based sampling method will be investigated in the future. For the following work, additional fission systems will be investigated with the proposed uncertainty propagation scheme. And a sensitivity analysis of effective multiplication factor to fission yields should be conducted to determine the reason behind the formation of effective multiplication factor skewed distribution.



---

*Source: 8014521-2020-09-25.xml*