Bayesian statistical model for cladding high-temperature burst under loss-of-coolant accident conditions

ABSTRACT

This study developed a probabilistic determination model with respect to cladding high-temperature burst conditions based on the Bayesian statistical method to reasonably evaluate fuel behaviors under loss-of-coolant accident conditions, including fuel fragmentation, relocation, and dispersal. The candidate models (correlations) were based on the widely accepted empirical model established based on nonirradiated fuel cladding data. Explanatory variables were added to improve the applicability of these models with respect to irradiated materials and generalization performance. The posterior predictive distribution of each candidate model was evaluated using Bayesian estimation comprising 238 sets of high-temperature burst test data. The generalization performance was evaluated using information criteria. The results of model evaluation showed improved predictive performance by considering the effect of hydrogen content. A comparison with burnup as an alternative explanatory variable confirmed that hydrogen content was the better parameter and other burnup-associated effects, such as irradiation hardening of the metal matrix and oxide growth (reduction of the metal matrix), were less dominant under burst conditions.

KEYWORDS:

• LOCA
• Zircaloy
• uncertainty

1. Introduction

When a loss-of-coolant accident (LOCA) occurs owing to a breach in the pressure boundary in a light water reactor, some fuel rods may swell and burst and the cladding tube forms an opening because of rapid reduction in coolant mass and pressure. This reduction results in the increase in fuel, cladding temperature, and rod internal pressure and decrease in cladding strength. In addition, for high-burnup fuel, a phenomenon known as fuel fragmentation, relocation, and dispersal (FFRD) has been observed [1,2]. Depending on conditions, such as burnup, irradiation history, and the relationship between the opening and fragmented pellet sizes, fragmented pellets may be released outside the fuel rods and/or deposited inside the cladding ballooning part [3]. The released fragmented pellets may be deposited on the grid spacer or lower structure of the fuel assembly. Also, the fuel temperature of the ballooning part filled with fragmented pellets can be higher [4]. Overall, FFRD may affect the coolability of the reactor core even under accident conditions that do not undergo guillotine breaks owing to embrittlement. Therefore, the importance of evaluation methodology with respect to the amount of released pellets in determining the impact of FFRD on the coolability is widely accepted and studied [5].

To evaluate the amount of fragmented pellets released outside the fuel rods, various factors related to FFRD, including cladding burst occurrence, cladding opening size, cladding deformation, pellet fragmentation occurrence, fragmented pellet size, pellet restraint by cladding (suppression of fragmentation), and relocation of fragmented pellets, etc., must be studied and modeled. However, the limitations related to experimental data and the large scatter observed in the existing data make assigning an accurate model for each factor challenging.

Conservative models are frequently employed to address these data limitations and uncertainties [6,7]. However, using these models can result in excessive conservatism, specifically when the evaluation involves multiple model components such as FFRD. Furthermore, the output of overly conservative estimations based on over-simplified models is typically less informative regarding physics/mechanisms underlying increased risk and, thus, provides less effective feedback with respect to fuel design and plant operations.

Another approach for constructing a best-estimate model, reported in a previous study [7,8], to fit the existing data, typically with maximum likelihood estimation, may be an alternative to simple conservative models. The maximum likelihood estimation method is a kind of point estimations, which originally does not consider uncertainties. Some distribution should be assumed to be necessary for an interval estimation; however, it is hard to approach the distribution in some cases. For example, if sufficient experimental data exist, uncertainty distribution involved in the experimental data can be considered as Student’s t-distribution. However, capturing the uncertainty trend based on a limited experimental dataset is still challenging, as in the case of FFRD and related behaviors. Additionally, the result obtained based on maximum likelihood estimation heavily relies on the corresponding fitted data. Consideration of generalization performance and the accuracy of the fitted data are often overlooked.

In modeling fuel behavior factors related to FFRD, avoiding overfitting (avoiding using many coefficients and variables) and naturally incorporating the effects of data sparseness and scattering are necessary. Therefore, uncertainty distribution in the data should not be definitive; for example, the parameters determining the uncertainty distribution can also be described as distributions, which allows the regression equation to have distribution of uncertainty. The component model must also be applicable to code-to-code coupling operations for a core-scale calculation. The integrated system of component codes and models can then realize a predictive distribution of the whole amount of released pellets from a reactor core to determine coolability with reasonably determined uncertainty having high robustness. The Bayesian statistical method is a practical method for developing component models, whose uncertainty distribution is automatically obtained using model output.

This study develops a probabilistic evaluation model for predicting cladding burst conditions as a crucial component model needed for integrated FFRD evaluations based on the Bayesian estimation. Some candidate models were developed, and the best model was chosen based on certain factors. One candidate model was created using a well-known empirical burst model established based on nonirradiated fuel cladding data. The other models were also created based on the empirical model and aimed to consider the effects of hydrogen content and burnup as additional explanatory variables. The model parameters were estimated using Bayesian estimation, supported by numerous data points obtained based on high-temperature cladding burst tests. The performances of the developed models were evaluated using information criteria to assess their generalization performance and overfitting.

2. Method

2.1 Candidate models

Numerous studies related to the burst behaviors of fuel cladding in LOCA conditions have been published [9]. High-temperature burst conditions in claddings are effectively described using functions of engineering hoop stress related to cladding at burst and cladding surface temperature at burst (Figure 1).

Figure 1. Burst temperature versus engineering hoop stress at burst related to cladding burst tests performed in studies [8] and [10–20] and the ORNL correlation.

The blue plots in Figure 1 show 238 sets of high-temperature burst tests. The effects of ramp rates, heating conditions, and atmospheres have been examined in preceding studies; however, data regarding high-burnup materials is insufficient [8] and [10–20]. Chapman of Oak Ridge National Laboratory (ORNL) [8] developed a burst temperature correlation derived from test data related to nonirradiated materials, which is also a function of ramp rate, as expressed in (Eq. 1(1) $\mathit{T}=3960-\frac{20.4\mathit{\varsigma }}{1+\frac{\mathit{R}}{28}}-\frac{8510000\mathit{\varsigma }}{100\left(1+\frac{\mathit{R}}{28}\right)+2790\mathit{\varsigma }}\left(0\le \mathit{R}\le 28\right),$(1) )

(1) $\mathit{T}=3960-\frac{20.4\mathit{\varsigma }}{1+\frac{\mathit{R}}{28}}-\frac{8510000\mathit{\varsigma }}{100\left(1+\frac{\mathit{R}}{28}\right)+2790\mathit{\varsigma }}\left(0\le \mathit{R}\le 28\right),$(1)

where $\mathit{T}\mathit{\hspace{0.17em}}\left({ }^{\circ }\mathrm{C}\right)$ is the surface temperature of cladding at burst, $\mathit{\varsigma }\mathit{\hspace{0.17em}}\left(\mathrm{kpsi}\right)$ is the engineering hoop stress of cladding at burst, and $\mathit{R}\mathit{\hspace{0.17em}}\left({ }^{\circ }\mathrm{C}/\mathrm{s}\right)$ is the ramp rate, with an applicable range of 0–28. The three lines shown in Figure 1 represent the model prediction of Eq. 1(1) $\mathit{T}=3960-\frac{20.4\mathit{\varsigma }}{1+\frac{\mathit{R}}{28}}-\frac{8510000\mathit{\varsigma }}{100\left(1+\frac{\mathit{R}}{28}\right)+2790\mathit{\varsigma }}\left(0\le \mathit{R}\le 28\right),$(1) with respect to the ramp rates of 0, 10, and 28.

While this model exhibits a specific performance in estimating the burst condition, the remaining discrepancy between the experimental data points and correlations can be reduced further. For instance, studies conducted at the Argonne National Laboratory (ANL) indicated that a high level of hydrogen content in the fuel cladding tubes of Zry-2 and ZIRLO (one of advanced Zr-based alloy fuel cladding tubes) affects burst conditions [21,22]. These studies by ANL include two kinds of non-irradiated samples: as-received and prehydrided fuel cladding tubes. The burst tests in their studies were thus conducted with as-received, prehydrided and irradiated fuel cladding tubes. The Japan Atomic Energy Agency (JAEA) studies confirm a similar trend using Zry-2, Zry-4 and some kinds of advanced Zr-based alloy fuel cladding tubes [23–25] (Figure 2). The pressure information in Figure 2 means the internal rod pressure at the beginning of each burst test. Figure 2 shows that hydrogen content affects the burst temperature. The remaining discrepancy observed in Figure 1 can be considerably reduced by introducing hydrogen content as an explanatory variable in Bayesian estimation, especially if sufficient supporting data becomes available. The abovementioned studies do not provide direct input to this study, which used cladding hoop stress level at the burst occurrence as the primary figure of merit.

Figure 2. Burst temperature versus hydrogen content in cladding burst tests performed by ANL and JAEA [21–25].

(1) The origin-based ORNL correlation

This study used Bayesian estimation in cases of some candidate models to predict burst conditions prepared based on Eq. 1(1) $\mathit{T}=3960-\frac{20.4\mathit{\varsigma }}{1+\frac{\mathit{R}}{28}}-\frac{8510000\mathit{\varsigma }}{100\left(1+\frac{\mathit{R}}{28}\right)+2790\mathit{\varsigma }}\left(0\le \mathit{R}\le 28\right),$(1) , with additional variables. Because using the mechanical response of the cladding tube as the output of the prediction model is more straightforward than using the burst temperature in the context of LOCA analyses, Eq. 1(1) $\mathit{T}=3960-\frac{20.4\mathit{\varsigma }}{1+\frac{\mathit{R}}{28}}-\frac{8510000\mathit{\varsigma }}{100\left(1+\frac{\mathit{R}}{28}\right)+2790\mathit{\varsigma }}\left(0\le \mathit{R}\le 28\right),$(1) is first solved for the cladding hoop stress at burst, with a slight change in the formulation to realize numerical stability. The unit of the engineering hoop stress of the cladding at burst, $\mathit{\varsigma }\left(\mathrm{kpsi}\right)$, is converted to σ (MPa) for convenience, and the coefficients and exponents of the factors (T and R, respectively) are partially simplified, as described in Eq. 2(2) $\mathit{\sigma }=C_1\times R^{0.5}+C_2\times T^{0.5}+C_3\times \mathit{T}+\mathit{I},$(2) .

(2) $\mathit{\sigma }=C_1\times R^{0.5}+C_2\times T^{0.5}+C_3\times \mathit{T}+\mathit{I},$(2)

where the coefficients C_1–3 and the intercept I are to be estimated using the Bayesian estimation. The intercept in this Bayesian statistical process effectively works as a residual error term, which means that it is challenging for the error to be reduced by other terms, including T and R because the error is not caused by a physics-based reason, for example. Eq. 2(2) $\mathit{\sigma }=C_1\times R^{0.5}+C_2\times T^{0.5}+C_3\times \mathit{T}+\mathit{I},$(2) is further refined to Eq. 3(3) $\mathit{\sigma }=C_1{\left(\mathit{R}-{\mathit{R}}_{\mathit{min}}\right)}^{0.5}+C_2{\left(\mathit{T}-{\mathit{T}}_{\mathit{min}}\right)}^{0.5}+C_3\left(\mathit{T}-T_{\mathit{min}}\right)+\mathit{I},$(3) (the origin-based ORNL correlation), which is one of the candidate models.

(3) $\mathit{\sigma }=C_1{\left(\mathit{R}-{\mathit{R}}_{\mathit{min}}\right)}^{0.5}+C_2{\left(\mathit{T}-{\mathit{T}}_{\mathit{min}}\right)}^{0.5}+C_3\left(\mathit{T}-T_{\mathit{min}}\right)+\mathit{I},$(3)

where all the explanatory variables in Eq. 3(3) $\mathit{\sigma }=C_1{\left(\mathit{R}-{\mathit{R}}_{\mathit{min}}\right)}^{0.5}+C_2{\left(\mathit{T}-{\mathit{T}}_{\mathit{min}}\right)}^{0.5}+C_3\left(\mathit{T}-T_{\mathit{min}}\right)+\mathit{I},$(3) (R and T) are calculated using the z-score, which is a scaled value that ensures the mean value of each variable is zero and the standard deviation (SD) is one for computational stability. Then, depending on the data, the values of the variables can be negative. The minimum values of ramp rate $R_{\mathit{min}}\mathit{\hspace{0.17em}}\left({ }^{\circ }\mathrm{C}/\mathrm{s}\right)$, and burst temperature $T_{\mathit{min}}\mathit{\hspace{0.17em}}\left({ }^{\circ }\mathrm{C}\right)$ are introduced to avoid the complicated calculations related to the exponent of negative values with decimals.

(2) The modified ORNL correlation combined with the FRAPCON-3.4 model

We consider other candidate models to develop a better model with reduced uncertainty and avoid overfitting. Since Eq. 3(3) $\mathit{\sigma }=C_1{\left(\mathit{R}-{\mathit{R}}_{\mathit{min}}\right)}^{0.5}+C_2{\left(\mathit{T}-{\mathit{T}}_{\mathit{min}}\right)}^{0.5}+C_3\left(\mathit{T}-T_{\mathit{min}}\right)+\mathit{I},$(3) refers to the equation derived initially from nonirradiated material and does not cover the effect of hydrogen content, it is natural for the newly developed model to increase the number of explanatory variables. Therefore, second candidate model includes hydrogen content $\mathit{H}\left(\mathrm{weight}\mathrm{parts}\mathrm{per}\left(\mathrm{wppm}\right)\right)$ in addition to Eq. 3(3) $\mathit{\sigma }=C_1{\left(\mathit{R}-{\mathit{R}}_{\mathit{min}}\right)}^{0.5}+C_2{\left(\mathit{T}-{\mathit{T}}_{\mathit{min}}\right)}^{0.5}+C_3\left(\mathit{T}-T_{\mathit{min}}\right)+\mathit{I},$(3) , called the modified ORNL correlation (Eq. 4(4) $\mathit{\sigma }=C_4{\left(\mathit{R}-{\mathit{R}}_{\mathit{min}}\right)}^{0.5}+C_5{\left(\mathit{T}-{\mathit{T}}_{\mathit{min}}\right)}^{0.5}+C_6\left(\mathit{T}-T_{\mathit{min}}\right)+\mathit{I}+C_7\mathit{H}.$(4) ) in this study,

(4) $\mathit{\sigma }=C_4{\left(\mathit{R}-{\mathit{R}}_{\mathit{min}}\right)}^{0.5}+C_5{\left(\mathit{T}-{\mathit{T}}_{\mathit{min}}\right)}^{0.5}+C_6\left(\mathit{T}-T_{\mathit{min}}\right)+\mathit{I}+C_7\mathit{H}.$(4)

Since the database compiled in this study includes only a small portion of datasets with H, we apply some well-established empirical correlations describing burnup-dependent hydrogen absorption. The high-temperature burst tests shown in Figure 1 primarily comprise Zircaloy-2 and −4 (Zry-2 and −4) cladding tubes. Although hydrogen absorption for Zry-4 is generally described as the combination of the oxide layer growth and hydrogen pick-up ratio, such a two-step approach adds complexity and, thus, the convergence problem in Bayesian estimation process. The simple FRAPCON-3.4 model developed for Zry-2 before 1998 (Eq. 5(5) $\begin{array}{rl}& \mathit{H}=47.8\exp \left(-\frac{1.3}{1+\mathit{G}}\right)+0.316\mathit{G}\mathit{where}\mathit{G}<50,\\ & \mathit{H}=28.9+\exp \left(0.117\left(\mathit{G}-20\right)\right)\mathit{where}\mathit{G}>50.\end{array}$(5) ) [26], which is a function of burnup $\mathit{G}$ (GWd/tU), is hence applied,

(5) $\begin{array}{rl}& \mathit{H}=47.8\exp \left(-\frac{1.3}{1+\mathit{G}}\right)+0.316\mathit{G}\mathit{where}\mathit{G}<50,\\ & \mathit{H}=28.9+\exp \left(0.117\left(\mathit{G}-20\right)\right)\mathit{where}\mathit{G}>50.\end{array}$(5)

Figure 3 compares Eq. 5(5) $\begin{array}{rl}& \mathit{H}=47.8\exp \left(-\frac{1.3}{1+\mathit{G}}\right)+0.316\mathit{G}\mathit{where}\mathit{G}<50,\\ & \mathit{H}=28.9+\exp \left(0.117\left(\mathit{G}-20\right)\right)\mathit{where}\mathit{G}>50.\end{array}$(5) and the data given in Figure 2, for which the prehydrided sample data are excluded from the comparison. Eq. 5(5) $\begin{array}{rl}& \mathit{H}=47.8\exp \left(-\frac{1.3}{1+\mathit{G}}\right)+0.316\mathit{G}\mathit{where}\mathit{G}<50,\\ & \mathit{H}=28.9+\exp \left(0.117\left(\mathit{G}-20\right)\right)\mathit{where}\mathit{G}>50.\end{array}$(5) roughly follows the experimental data trends, including Zry-2 and −4 data. The burnup $\mathit{G}$ was assumed to be zero for a few datasets whose sample preparation conditions were not identified from the literature.

Figure 3. Comparison between the FRAPCON-3.4 model and as-received and irradiated cladding tubes.

For a linear regression with parametric Bayesian estimation, this study assumes that the error ϵ between the solution of the regression equation σ and the measured value ${\sigma }_{\mathit{obs}}$ follows a normal distribution. The relationship is then expressed in Eqs. 6(6) ${\sigma }_{\mathit{obs}}=\mathit{\sigma }+\mathit{ϵ},$(6) and 7(7) $\mathit{ϵ}\sim \mathit{Normal}\left(0,\mathit{SD}\right)$(7) ,

(6) ${\sigma }_{\mathit{obs}}=\mathit{\sigma }+\mathit{ϵ},$(6)

(7) $\mathit{ϵ}\sim \mathit{Normal}\left(0,\mathit{SD}\right)$(7)

where $\mathit{SD}$ is standard deviation. Namely, the distribution of the experimental results follows as Eq. 8(8) ${\sigma }_{\mathit{obs}}\sim \mathit{Normal}\left(\mathit{\sigma },\mathit{SD}\right).$(8) ,

(8) ${\sigma }_{\mathit{obs}}\sim \mathit{Normal}\left(\mathit{\sigma },\mathit{SD}\right).$(8)

(3) The nonparametric model with the ramp rate, burst temperature, and hydrogen content

The nonparametric model defined as Eq. 9(9) $\mathit{\sigma }=\mathit{GP}\left(\mathit{R},\mathit{T},\mathit{H}\right).$(9) does not have a specific equation form like Eq. 3(3) $\mathit{\sigma }=C_1{\left(\mathit{R}-{\mathit{R}}_{\mathit{min}}\right)}^{0.5}+C_2{\left(\mathit{T}-{\mathit{T}}_{\mathit{min}}\right)}^{0.5}+C_3\left(\mathit{T}-T_{\mathit{min}}\right)+\mathit{I},$(3) . This nonparametric model was used to verify the possibility that the formulation of Eq. 3(3) $\mathit{\sigma }=C_1{\left(\mathit{R}-{\mathit{R}}_{\mathit{min}}\right)}^{0.5}+C_2{\left(\mathit{T}-{\mathit{T}}_{\mathit{min}}\right)}^{0.5}+C_3\left(\mathit{T}-T_{\mathit{min}}\right)+\mathit{I},$(3) might already be overfitting for nonirradiated material and deteriorate the prediction performance for irradiated material. The number of explanatory variables in Eq. 9(9) $\mathit{\sigma }=\mathit{GP}\left(\mathit{R},\mathit{T},\mathit{H}\right).$(9) is the same as in Eq. 4(4) $\mathit{\sigma }=C_4{\left(\mathit{R}-{\mathit{R}}_{\mathit{min}}\right)}^{0.5}+C_5{\left(\mathit{T}-{\mathit{T}}_{\mathit{min}}\right)}^{0.5}+C_6\left(\mathit{T}-T_{\mathit{min}}\right)+\mathit{I}+C_7\mathit{H}.$(4) , and the Gaussian process, described later, performs the estimation,

(9) $\mathit{\sigma }=\mathit{GP}\left(\mathit{R},\mathit{T},\mathit{H}\right).$(9)

The above formulas of Eqs. 3(3) $\mathit{\sigma }=C_1{\left(\mathit{R}-{\mathit{R}}_{\mathit{min}}\right)}^{0.5}+C_2{\left(\mathit{T}-{\mathit{T}}_{\mathit{min}}\right)}^{0.5}+C_3\left(\mathit{T}-T_{\mathit{min}}\right)+\mathit{I},$(3) and 4(4) $\mathit{\sigma }=C_4{\left(\mathit{R}-{\mathit{R}}_{\mathit{min}}\right)}^{0.5}+C_5{\left(\mathit{T}-{\mathit{T}}_{\mathit{min}}\right)}^{0.5}+C_6\left(\mathit{T}-T_{\mathit{min}}\right)+\mathit{I}+C_7\mathit{H}.$(4) are the basic candidate models, and Eq. 9(9) $\mathit{\sigma }=\mathit{GP}\left(\mathit{R},\mathit{T},\mathit{H}\right).$(9) is the additional verification model for a critical assessment.

The nonparametric model (Eq. 9(9) $\mathit{\sigma }=\mathit{GP}\left(\mathit{R},\mathit{T},\mathit{H}\right).$(9) ) adopts the Gaussian process using a multivariate normal distribution (Eq. 10(10) $\mathit{P}\left(\mathit{X}\right)\sim \mathit{MultiNormal}\left(0,\mathit{\kappa }\left(\mathit{X},{\mathit{X}}^{\prime }\right)\right),$(10) ). The factor X in Eq. 10(10) $\mathit{P}\left(\mathit{X}\right)\sim \mathit{MultiNormal}\left(0,\mathit{\kappa }\left(\mathit{X},{\mathit{X}}^{\prime }\right)\right),$(10) represents the input vector at which a prediction is made, and X′ represents the input vectors from the training dataset. The Gaussian kernel formula $\mathit{\kappa }\left(\mathit{X},{\mathit{X}}^{\prime }\right)$ shown in Eq. 11(11) $\mathit{\kappa }\left(\mathit{X},{\mathit{X}}^{{ }^{\mathbf{\prime }}}\right)={\theta }_{1}exp\left(-\frac{{\left(\mathit{X}-{\mathit{X}}^{\prime }\right)}^2}{2{\theta }_2}\right)+\mathit{diag}\left({\theta }_3\right).$(11) describes the similarity between the data ($\mathit{X}$ and ${\mathit{X}}^{\prime }$), which is the essence of the Gaussian process to give an infinite-dimensional Gaussian distribution for the input. The parameters in the kernel function θ_1,2,3 are also to be calculated in the Bayesian estimation.

(10) $\mathit{P}\left(\mathit{X}\right)\sim \mathit{MultiNormal}\left(0,\mathit{\kappa }\left(\mathit{X},{\mathit{X}}^{\prime }\right)\right),$(10)

(11) $\mathit{\kappa }\left(\mathit{X},{\mathit{X}}^{{ }^{\mathbf{\prime }}}\right)={\theta }_{1}exp\left(-\frac{{\left(\mathit{X}-{\mathit{X}}^{\prime }\right)}^2}{2{\theta }_2}\right)+\mathit{diag}\left({\theta }_3\right).$(11)

The Bayesian estimation method first requires prior distributions to each parameter, and the distributions are updated based on the database (Bayesian update) to obtain the posterior distributions. The parameters to be estimated using the Bayesian method are $C_{1-7}$ in Eqs. 3(3) $\mathit{\sigma }=C_1{\left(\mathit{R}-{\mathit{R}}_{\mathit{min}}\right)}^{0.5}+C_2{\left(\mathit{T}-{\mathit{T}}_{\mathit{min}}\right)}^{0.5}+C_3\left(\mathit{T}-T_{\mathit{min}}\right)+\mathit{I},$(3) and 4(4) $\mathit{\sigma }=C_4{\left(\mathit{R}-{\mathit{R}}_{\mathit{min}}\right)}^{0.5}+C_5{\left(\mathit{T}-{\mathit{T}}_{\mathit{min}}\right)}^{0.5}+C_6\left(\mathit{T}-T_{\mathit{min}}\right)+\mathit{I}+C_7\mathit{H}.$(4) , I, SD, and ${\theta }_{1,2,3}$ in Eq. 11(11) $\mathit{\kappa }\left(\mathit{X},{\mathit{X}}^{{ }^{\mathbf{\prime }}}\right)={\theta }_{1}exp\left(-\frac{{\left(\mathit{X}-{\mathit{X}}^{\prime }\right)}^2}{2{\theta }_2}\right)+\mathit{diag}\left({\theta }_3\right).$(11) . This study assumes all prior distributions of $C_{1-7}$ in Eqs. 3(3) $\mathit{\sigma }=C_1{\left(\mathit{R}-{\mathit{R}}_{\mathit{min}}\right)}^{0.5}+C_2{\left(\mathit{T}-{\mathit{T}}_{\mathit{min}}\right)}^{0.5}+C_3\left(\mathit{T}-T_{\mathit{min}}\right)+\mathit{I},$(3) and 4(4) $\mathit{\sigma }=C_4{\left(\mathit{R}-{\mathit{R}}_{\mathit{min}}\right)}^{0.5}+C_5{\left(\mathit{T}-{\mathit{T}}_{\mathit{min}}\right)}^{0.5}+C_6\left(\mathit{T}-T_{\mathit{min}}\right)+\mathit{I}+C_7\mathit{H}.$(4) , I, SD, and ${\theta }_{1,2,3}$ in Eq. 11(11) $\mathit{\kappa }\left(\mathit{X},{\mathit{X}}^{{ }^{\mathbf{\prime }}}\right)={\theta }_{1}exp\left(-\frac{{\left(\mathit{X}-{\mathit{X}}^{\prime }\right)}^2}{2{\theta }_2}\right)+\mathit{diag}\left({\theta }_3\right).$(11) to be an uninformative uniform distribution. Since many kinds of experimental errors generally follow Gaussian distributions, using the distributions as the prior distributions is likely valid and additionally has the advantage that iterative calculations during Bayesian estimation are more likely to converge. Meanwhile, the uniform distribution is a more careful or conservative approach because it is the assumption that best eliminates the subjectivity and bias of an estimator. The uniform distribution is chosen because it is the first time to apply the Bayesian estimation method to determine cladding burst condition, and convergence problems didn’t occur even using the uniform distribution.

Since it is generally difficult to obtain the analytical results of integral calculations associated with the Bayesian updating, the No-U-Turn Sampler (NUTS) method, a Markov Chain Monte Carlo (MCMC) method, becomes a solution [27]. The specific computational tool is the programming language Stan, which implements the NUTS method as a standard algorithm and is called from PyStan, a library for Python [28]. The identical database was applied to the candidate models (Eqs. 3(3) $\mathit{\sigma }=C_1{\left(\mathit{R}-{\mathit{R}}_{\mathit{min}}\right)}^{0.5}+C_2{\left(\mathit{T}-{\mathit{T}}_{\mathit{min}}\right)}^{0.5}+C_3\left(\mathit{T}-T_{\mathit{min}}\right)+\mathit{I},$(3) and 4(4) $\mathit{\sigma }=C_4{\left(\mathit{R}-{\mathit{R}}_{\mathit{min}}\right)}^{0.5}+C_5{\left(\mathit{T}-{\mathit{T}}_{\mathit{min}}\right)}^{0.5}+C_6\left(\mathit{T}-T_{\mathit{min}}\right)+\mathit{I}+C_7\mathit{H}.$(4) ) and the verification model (Eq. 9(9) $\mathit{\sigma }=\mathit{GP}\left(\mathit{R},\mathit{T},\mathit{H}\right).$(9) ).

2.2 Evaluation method

The widely applicable information criterion (WAIC) [29] and the leave-one-out cross-validation (LOOCV) [30] methods are adopted to evaluate the generalization performance of the candidate models. Although the widely applicable Bayesian information criterion (WBIC) is as popular as WAIC, WBIC is generally accepted as a measure of the fitness of a created model with respect to existing data. Using WBIC for evaluating the candidate models presented in this study is not desirable because most of the tests are performed using low burnup rods. Because the primary focus of this study is regarding the applicability of the models toward FFRD-related evaluation, evaluating the generalization performance with respect to the sparse data region at high burnup is necessary. The WAIC and LOOCV are thus better suited.

It is defined in the WAIC that a difference between a generalization loss (an inverse of the generalization performance) and an empirical loss E (Eq. 12(12) $\mathit{E}=-\frac{1}{n}\sum _{\mathit{i}=1}^{n}log{E}_w\left[\mathit{p}\left(X_i|\mathit{w}\right)\right],$(12) ) can be expressed as the functional variance V (Eq. 13(13) $\mathit{V}=\sum _{\mathit{i}=1}^n\left\{E_w\left[{\left(logp\left(X_i|\mathit{w}\right)\right)}^2\right]-E_w{\left[logp\left(X_i|\mathit{w}\right)\right]}^2\right\},$(13) ) divided by the number of data n. The WAIC is shown in Eq. 14(14) $\mathit{WAIC}=\left(\mathit{E}+\frac{V}{n}\right),$(14) ,

(12) $\mathit{E}=-\frac{1}{n}\sum _{\mathit{i}=1}^{n}log{E}_w\left[\mathit{p}\left(X_i|\mathit{w}\right)\right],$(12)

(13) $\mathit{V}=\sum _{\mathit{i}=1}^n\left\{E_w\left[{\left(logp\left(X_i|\mathit{w}\right)\right)}^2\right]-E_w{\left[logp\left(X_i|\mathit{w}\right)\right]}^2\right\},$(13)

(14) $\mathit{WAIC}=\left(\mathit{E}+\frac{V}{n}\right),$(14)

where $\mathit{p}$ is the predictive distribution, $X_i$ is the data, $\mathit{w}$ is the parameter subjected to the Bayesian estimation, and $E_w$ is the expected value over $\mathit{w}$. The empirical loss means the unfitness of the created model for the existing data, and the functional variance divided by n means the effect of uncertainties caused by unknown data.

In LOOCV, one test is once removed from the entire test, and all remaining tests are used for training. Namely, the omitted one is regarded as unknown and used for evaluating the fitness of the trained model. Each test receives the removal treatment, and the average value of such evaluations is defined as the LOOCV value. This study uses the importance-sampling cross-validation: an approximated LOOCV method frequently used in a Bayesian estimation for reducing computation time [31]. Eq. 15(15) $\mathit{LOOCV}=\frac{1}{n}\sum _{\mathit{i}=1}^{n}log{E}_w\left[1/\mathit{p}\left(X_i|\mathit{w}\right)\right].$(15) expresses the approximated method,

(15) $\mathit{LOOCV}=\frac{1}{n}\sum _{\mathit{i}=1}^{n}log{E}_w\left[1/\mathit{p}\left(X_i|\mathit{w}\right)\right].$(15)

We evaluated the WAIC and LOOCV for double-checking. Lower values indicate increased generalization performance. The WAIC and LOOCV evaluated values are known to become asymptotically equivalent [29]. If the two evaluations give far different results, there could be problems with the quality and/or quantity of data.

3. Results & discussion

The convergence of the posterior distributions estimated using MCMC sampling was confirmed based on the Gelman – Rubin statistic (Rhat) [32], which is also employed in PyStan. When Rhat is close to 1.0, it means that the Markov Chain properly converges to a steady state. Herein, all parameters in each model subjected to Bayesian estimation were equal to 1.0 at Rhat; thus, this study could realize valid posterior distributions.

3.1. Uncertainty distributions

Figs. 4–9 demonstrate the posterior predictive distribution of engineering hoop stress at burst for each model obtained based on the processes described in the previous section. Figs. 5 and 7 exhibit the two-dimensional distributions for the ramp rate 5°C/s, which is a typical condition related to a postulated LOCA sequence [33].

Figure 4. Three-dimensional distribution of uncertainty related to the origin-based ORNL correlation.

Figure 5. Two-dimensional distribution of uncertainty related to the origin-based correlation, wherein the fixed ramp rate is 5 (°C/s).

Figure 6. Three-dimensional distribution of uncertainty related to the modified ORNL correlation wherein the fixed hydrogen content is 100 (wppm).

Figs. 4, 6 and 9 show the transparent red, green, and blue surfaces, i.e. the 97.5th percentile, the median, and the 2.5th percentile of the predictive distribution of engineering hoop stress of cladding at burst, respectively. Namely, the area bounded by the red and blue surfaces is the 95% predictive interval: the uncertainty contained in each model. The blue points are the literature data.

The engineering hoop stress of cladding at burst obtained by Bayesian estimation has uncertainty distributions. It reflects the sparseness/density of the datasets depending on the conditions defined by the explanatory variables, providing robustness in uncertainty quantification and giving larger uncertainties for the range less supported by the data.

According to the predictive intervals in Figs. 5 and 7, the burst stress predicted by the modified ORNL correlation (Eq. 4(4) $\mathit{\sigma }=C_4{\left(\mathit{R}-{\mathit{R}}_{\mathit{min}}\right)}^{0.5}+C_5{\left(\mathit{T}-{\mathit{T}}_{\mathit{min}}\right)}^{0.5}+C_6\left(\mathit{T}-T_{\mathit{min}}\right)+\mathit{I}+C_7\mathit{H}.$(4) ) with hydrogen content of 100 (wppm) is lower than that predicted by the origin-based correlation (Eq. 3(3) $\mathit{\sigma }=C_1{\left(\mathit{R}-{\mathit{R}}_{\mathit{min}}\right)}^{0.5}+C_2{\left(\mathit{T}-{\mathit{T}}_{\mathit{min}}\right)}^{0.5}+C_3\left(\mathit{T}-T_{\mathit{min}}\right)+\mathit{I},$(3) ). It confirms that the effect of hydrogen content on the burst condition is reflected on the obtained Bayesian model as intended. The two-dimensional plot of cladding hoop stress at burst in Figure 7 visualizes the large uncertainty band more quantitatively: 10–70 MPa, for example, for the surface temperature (800°C).

Figure 7. Two-dimensional distribution of uncertainty related to the modified ORNL correlation wherein the fixed ramp rate is 5 (°C/s) and hydrogen content is 100 (wppm).

Figure 8 visualizes the sensitivity of the modified ORNL correlation (Eq. 4(4) $\mathit{\sigma }=C_4{\left(\mathit{R}-{\mathit{R}}_{\mathit{min}}\right)}^{0.5}+C_5{\left(\mathit{T}-{\mathit{T}}_{\mathit{min}}\right)}^{0.5}+C_6\left(\mathit{T}-T_{\mathit{min}}\right)+\mathit{I}+C_7\mathit{H}.$(4) ) to hydrogen content. The hydrogen content of 13 (wppm) corresponds to the minimum level given by Eq. 5(5) $\begin{array}{rl}& \mathit{H}=47.8\exp \left(-\frac{1.3}{1+\mathit{G}}\right)+0.316\mathit{G}\mathit{where}\mathit{G}<50,\\ & \mathit{H}=28.9+\exp \left(0.117\left(\mathit{G}-20\right)\right)\mathit{where}\mathit{G}>50.\end{array}$(5) . Fig. 8 indicates that the model predicts earlier bursts for higher hydrogen contents, which should be the primary explanatory variable to represent the burnup effect. For the hydrogen content higher than 200 (wppm), the 95% predictive interval estimated by the modified ORNL correlation is systematically lower than the original correlation value. Namely, with the modified ORNL correlation, a high-burnup state will be processed as more penalized regarding, for example, FFRD.

Figure 8. The change in uncertainty by fixing hydrogen content with respect to the modified ORNL correlation wherein the fixed ramp rate is 5 (°C/s).

As the hydrogen effect is incorporated in the modified ORNL correlation as a linear term (Eq. 4(4) $\mathit{\sigma }=C_4{\left(\mathit{R}-{\mathit{R}}_{\mathit{min}}\right)}^{0.5}+C_5{\left(\mathit{T}-{\mathit{T}}_{\mathit{min}}\right)}^{0.5}+C_6\left(\mathit{T}-T_{\mathit{min}}\right)+\mathit{I}+C_7\mathit{H}.$(4) ), this formulation entails the monotonic tendency of hydrogen effect in any data region. However, since there is no data directly supporting the model for high hydrogen content, the predictive intervals shown in Figure 8 get spread in higher hydrogen content thanks to the Bayesian approach. Namely, according to the derived modified ORNL model, it becomes increasingly uncertain whether earlier bursts occur in higher hydrogen content, and the risk of model extrapolation to the data scarce (or lacking) range is compensated.

The 95% predictive interval in the nonparametric model (Figure 9) looks deformed in the data-sparse region, specifically the low ramp rate and high-temperature region. The region well-supported by data, specifically the low-temperature region, exhibits a narrower uncertainty band, but the surface shape is more skewed than those in the parametric models (Figures 4 and 6) in the temperature and ramp rate dimensions. Due to the extreme flexibility of the Gaussian process and the distribution of data sparseness, the nonparametric model fails to capture the physics-based trend, like the monotonic decrease in burst stress with increasing burst temperature, associated with material softening, as far as the present Bayesian inference reached. Although it is not essential because an actual ramp rate in an LOCA sequence is several °C/s, a capacity exists to investigate the reason and validity that might improve the uncertainty band in the parametric models by fixing the coefficients and components of the ramp rate.

Figure 9. Three-dimensional distribution of uncertainty with respect to the nonparametric model, wherein the fixed hydrogen content is 13 (wppm).

3.2 Generalization performance

Table 1 lists the evaluated values of the WAIC and LOOCV for each model. The model described by Eq. 4(4) $\mathit{\sigma }=C_4{\left(\mathit{R}-{\mathit{R}}_{\mathit{min}}\right)}^{0.5}+C_5{\left(\mathit{T}-{\mathit{T}}_{\mathit{min}}\right)}^{0.5}+C_6\left(\mathit{T}-T_{\mathit{min}}\right)+\mathit{I}+C_7\mathit{H}.$(4) performs better than Eq. 3(3) $\mathit{\sigma }=C_1{\left(\mathit{R}-{\mathit{R}}_{\mathit{min}}\right)}^{0.5}+C_2{\left(\mathit{T}-{\mathit{T}}_{\mathit{min}}\right)}^{0.5}+C_3\left(\mathit{T}-T_{\mathit{min}}\right)+\mathit{I},$(3) . On the other hand, the verification model (Eq. 9(9) $\mathit{\sigma }=\mathit{GP}\left(\mathit{R},\mathit{T},\mathit{H}\right).$(9) ) exhibits the best performance.

Table 1. Evaluated WAIC and LOOCV.

Model		WAIC	LOOCV
Eq. 3	The origin-based ORNL correlation	0.830	0.830
Eq. 4	The modified ORNL correlation combined with FRAPCON-3.4 model	0.739	0.739
	The modified ORNL correlation directly applying burnup instead of hydrogen content	0.765	0.765
Eq. 9	The nonparametric model with ramp rate, burst temperature and hydrogen content	0.513	0.513

The evaluation of the generalization performance of each candidate model using the WAIC and LOOCV showed improved performance of the modified ORNL correlation with an additional explanatory variable of hydrogen content (Eq. 4(4) $\mathit{\sigma }=C_4{\left(\mathit{R}-{\mathit{R}}_{\mathit{min}}\right)}^{0.5}+C_5{\left(\mathit{T}-{\mathit{T}}_{\mathit{min}}\right)}^{0.5}+C_6\left(\mathit{T}-T_{\mathit{min}}\right)+\mathit{I}+C_7\mathit{H}.$(4) ) than the origin-based ORNL model (Eq. 3(3) $\mathit{\sigma }=C_1{\left(\mathit{R}-{\mathit{R}}_{\mathit{min}}\right)}^{0.5}+C_2{\left(\mathit{T}-{\mathit{T}}_{\mathit{min}}\right)}^{0.5}+C_3\left(\mathit{T}-T_{\mathit{min}}\right)+\mathit{I},$(3) ). Hydrogen content was accordingly confirmed as crucial for improving the reliability of cladding burst determination.

Adopting burnup $\mathit{G}$ instead of $\mathit{H}$ as the additional explanatory variable in Eq. 4(4) $\mathit{\sigma }=C_4{\left(\mathit{R}-{\mathit{R}}_{\mathit{min}}\right)}^{0.5}+C_5{\left(\mathit{T}-{\mathit{T}}_{\mathit{min}}\right)}^{0.5}+C_6\left(\mathit{T}-T_{\mathit{min}}\right)+\mathit{I}+C_7\mathit{H}.$(4) resulted in less improvement of the generalization performance, as observed in Table 1. Hydrogen effect in the current database is not likely a spurious correlation by other unnoticed parameters, confirming the primary importance of hydrogen effect among various burnup-related effects, such as irradiation damage in the metal matrix and the reduction in the metal matrix by oxide growth.

According to the information criteria, the nonparametric model (Eq. 9(9) $\mathit{\sigma }=\mathit{GP}\left(\mathit{R},\mathit{T},\mathit{H}\right).$(9) ) showed highly increased performance compared to the parametric models (Eqs. 3(3) $\mathit{\sigma }=C_1{\left(\mathit{R}-{\mathit{R}}_{\mathit{min}}\right)}^{0.5}+C_2{\left(\mathit{T}-{\mathit{T}}_{\mathit{min}}\right)}^{0.5}+C_3\left(\mathit{T}-T_{\mathit{min}}\right)+\mathit{I},$(3) and 4(4) $\mathit{\sigma }=C_4{\left(\mathit{R}-{\mathit{R}}_{\mathit{min}}\right)}^{0.5}+C_5{\left(\mathit{T}-{\mathit{T}}_{\mathit{min}}\right)}^{0.5}+C_6\left(\mathit{T}-T_{\mathit{min}}\right)+\mathit{I}+C_7\mathit{H}.$(4) ). However, the substantial extension of the uncertainty band observed in the data-sparse region in Figure 9 did not have a negative impact on the WAIC and LOOCV evaluations.

As highlighted in the previous section, for the current burst test database, with its specific data sparseness and/or distribution of data sparseness, the Gaussian process fails to capture the trend expected from the underlying physics related to cladding high-temperature burst and produces less practical model, with a substantial uncertainty band. Prior knowledge, such as the temperature effect on cladding mechanical properties, explicitly used in the ORNL correlations, is essential to establish a reasonable model as a component of the integral LOCA simulation system.

Nevertheless, the nonparametric model (Eq. 9(9) $\mathit{\sigma }=\mathit{GP}\left(\mathit{R},\mathit{T},\mathit{H}\right).$(9) ) is an informative and a helpful model to conduct the sensitivity study, when implemented in an integrated simulation. This model can predict the evaluation risk by reacting as an especially large uncertainty band, as expected from the model response observed in Figure 9. Additionally, substantial reduction in the WAIC and LOOCV values related to the nonparametric model indicates that the form of ORNL correlation is nearly optimized to the behaviors of nonirradiated cladding samples, with room for improvement with respect to the increase in balance between the generalization performances related to nonirradiated and irradiated cladding samples.

4. Conclusions

This study developed a probabilistic determination model with respect to cladding burst conditions, incorporating data uncertainties using the Bayesian statistical method. The Bayesian inference method produced two parametric models based on the widely used ORNL correlation and one nonparametric model based on the Gaussian process, providing burst stress prediction with uncertainty distribution. A comparison between the two parametric models with and without hydrogen content as an explanatory variable using the WAIC and LOOCV methods showed improved reliability with respect to cladding burst determination. The sensitivity study related to incorporating burnup instead of hydrogen content showed that hydrogen content is the more relevant parameter for evaluating cladding bursts and other burnup-associated effects are not dominant. In the case that hydrogen content is more than 200 (wppm), the difference of predictions between the modified ORNL correlation and the original ORNL correlation was prominent. Since this level of hydrogen content is realistic in high burnup fuels, considering the hydrogen effect is crucial for evaluating FFRD. Although the nonparametric model was evaluated as the best model based on the information criteria, its uncertainty distribution is highly skewed to capture the physics-based trend. Thus, the modified ORNL correlation incorporated with the hydrogen content variable is a recommended correlation based on Bayesian models for developing a probabilistic burst model as a component of integral LOCA simulations, which is designed in numerous recent studies for evaluating FFRD. The burst model is well-supported by the experimental data for the temperature range of 700°C–900°C at a low ramp rate (~5°C/s). Further robustness with respect to the treatment of model uncertainty can be increased via the combined use of nonparametric and parametric models, which realize substantially large uncertainty responses in data-sparse conditions.

Acknowledgments

Part of the analyses in this study was performed under the research entrusted by the Nuclear Regulation Authority, Japan. The authors express their sincere gratitude to Mr. K. Hattori for performing data analysis based on the Bayesian inference method.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

MCMC samples of model parameters posterior distributions and the rotatable 3D graphs of Figs. 4, 6 and 9 are available on Mendeley online storage: https://data.mendeley.com/datasets/cm59yssdmb/2. The distributions of uncertainty shown in this study can be calculated using this dataset.