Uncertainty and Sensitivity Analysis of a Dry Cask for Spent Nuclear Fuel

Abstract

Nuclear safety relies to a good extent on thoroughly validated codes. However, code predictions are affected by uncertainties that need to be quantified for a more accurate evaluation of safety margins. In this regard, the present paper proposes a preliminary uncertainty and sensitivity analysis of the thermal behavior of a concrete-based dry cask for spent nuclear fuel storage, employing the MELCOR code and a series of MATLAB scripts. As thermal behavior is of utmost importance for the fulfillment of United States Nuclear Regulatory Commission (USNRC) safety requirements, the Peak Cladding Temperature (PCT) has been addressed as the key Figure of Merit (FOM). Variables related to the main heat transfer mechanisms have been selected as input parameters for the uncertainty quantification, whereas heat source and heat sink, namely decay power and external air temperature, have been dealt with in a separate sensitivity analysis. The results show that the selected parameters have a weak influence on the PCT, whereas it is strongly related to the decay power and external air temperature values. In any case, PCT stays below the regulatory threshold even under the considered off-normal conditions.

1. Introduction

The management of spent nuclear fuel (SNF) is recognized as the focal point of the back end of the nuclear fuel cycle. Both open and closed fuel cycle options have been taken into consideration by different countries, even though most of them have shown a greater tendency towards the former. However, delays in the introduction of permanent geological disposal facilities and in governments’ decisions about reprocessing processes have led to the need for storage systems for SNF other than spent fuel pools, which are reaching their maximum capacity. In this framework, dry cask systems have been proposed and accepted worldwide as an interim storage system for SNF for up to 100 years [1,2,3,4,5].

As other nuclear systems, dry casks have to fulfill the safety requirements as stated in United States Nuclear Regulatory Commission (USNRC) regulations [6]. In particular, to guarantee fuel rods’ integrity and because most degradation mechanisms are temperature-dependent, a threshold temperature of 673 K has been set [7] for fuel cladding under normal conditions (843 K for off-normal and accidental conditions). The above mentioned temperatures are conservatively derived in order to avoid creep and to limit the hoop stress in the cladding (the hoop stress value limit is set to 90 MPa to preclude hydride reorientation [7]). In this regard, the thermal analysis of the SNF has been identified as high-priority by US Department of Energy (DOE) [8,9].

In the past, many simulations have been carried out to study the dry cask behavior, mostly by employing thermal-hydraulic codes based on the finite volume method [10,11]. On the contrary, very few experimental tests have been performed: one experimental campaign was conducted by EPRI in 1987 for pressurized water reactor (PWR) spent fuel assemblies [11], while more recently a dry cask simulator (DCS) was built in Sandia National Laboratories (SNL) to reproduce boiling water reactor (BWR) spent fuel assemblies [12]. Moreover, heat removal tests were conducted in Japan on two types of full-scale concrete casks, taking into account different decay heat values [13]. Data from these experiments have been employed in several works [14,15,16,17,18,19].

In recent years, a great effort has been devoted to the thermal analysis of dry casks with computational fluid dynamic (CFD) codes [14,20,21,22,23,24,25]. However, despite the ability of CFD codes to provide deep insights, their high computational cost represents an important drawback that opens the way to simpler approaches [26]. Moreover, the computational burden would be an obstacle to the comprehensive implementation of uncertainty quantification (UQ) analyses. In fact, even though increasing attention has been paid to UQ in recent times [27,28], almost no proper uncertainty analysis (UA) on dry casks has been found in the literature. One example can be found in [29], where validation and uncertainty quantification for a CFD model of a metallic cask are proposed. Numerical and input uncertainties were investigated, and, as a result, an uncertainty of around 60 K was noted in the peak cladding temperature (PCT).

In this framework, the present paper reports the results of a preliminary uncertainty and sensitivity analysis (UASA) based on a MELCOR model of the HI-STORM 100S dry cask [25,30]. The analysis focuses on the PCT, being the variable under regulatory surveillance, and on the uncertainties embedded in the values predicted by the code. A reasonable accuracy and a significant reduction of the computational cost make the MELCOR code [31,32] suitable for performing UQ in an easier and faster way.

2. Modeling

2.1. System and Scenario Description

The HI-STORM (Holtec International Storage and Transport Operation Reinforced Model) system is a dry storage for SNF designed to provide confinement, radiation shielding, criticality control and heat removal in an autonomous way. It consists of a sealed metallic Multi-Purpose Canister (MPC)—backfilled with high-purity helium and vertically loaded with SNF—placed inside a concrete overpack. The system relies on the thermosiphon effect of the helium inside the canister combined with the natural circulation of the external air to passively cool the spent fuel, as shown in Figure 1.

The cask design addressed in this work is the HI-STORM 100S with a MPC-32 [25]. The metallic canister is a confinement system which consists of a stainless steel enclosure vessel equipped with a honeycombed fuel basket where 32 PWR 17 × 17 fuel assemblies are located. A dry inert atmosphere inside the MPC is guaranteed by the presence of pure helium as the backfill gas.

The overpack consists instead of a steel-layered concrete structure which provides both structural strength and radiological protection. Four circumferentially distributed inlet ducts at the bottom of the overpack allow the external air to enter the air gap between the canister and the overpack and to cool the MPC. Four outlet ducts at the top of the overpack allow the air to discharge the heat into the environment and to ensure effective natural circulation.

The main dimensions of the above-mentioned sub-systems are reported in

Table 1. The overpack main dimensions [25].

-	Dimension
Total Height (m)	5.824
Cavity Height (m)	5.22
Lid Thickness (m)	0.48
Baseplate Thickness (m)	0.124
MPC Support Thickness (m)	0.076
External Diameter (m)	3.37
Internal Diameter (m)	1.87
MPC Support Diameter (m)	1.65
Inlet Dimensions (m)	0.76 × 0.0762
Outlet Dimensions (m)	0.6477 × 0.12065
Total Inlet Area (m2)	0.232
Total Outlet Area (m2)	0.313

,

Table 2. The MPC-32 main dimensions [25].

-	Dimension
Total Height (m)	4.85
Cavity Height (m)	4.50
Lid Thickness (m)	0.29
Baseplate Thickness (m)	0.06
External Diameter (m)	1.78
Internal Diameter (m)	1.754
Cylindrical Wall Thickness (m)	0.013
Upper Fuel Spacer Height (m)	0.21
Lower Fuel Spacer Height (m)	0.20

and

Table 3. The honeycombed fuel basket main dimensions [25].

-	Dimension
Basket Height (m)	4.44
Cell Dimension (m)	0.22 × 0.22
Thickness (m)	0.0071
Distance to the MPC lid (m)	0.06
External Diameter (m)	1.78
Internal Diameter (m)	1.754
Opening Radius (m)	0.0635
Opening Area (m2)	0.00633

, respectively.

In particular, for the steady-state scenario addressed in this work, an initial temperature of 294.26 K and a pressure of 0.334 MPa have been taken into account for the helium inside the canister, whereas the temperature and pressure of external air have been taken equal to 300 K and 0.101 MPa. Additionally, the total decay power coming from the fuel assemblies has been considered equal to 30.0 kW. The initial and boundary conditions have been set as in [22].

2.2. MELCOR Model

The model for the HI-STORM 100S, better described in [30], is briefly summarized in this section. The MELCOR 2.2 code, which is a fully integrated computer code to model accident progression in LWRs, has been employed. Despite the fact that dry casks are completely outside MELCOR’s domain of application, the code’s structure in form of packages and its integral nature have been the key to model the entire HI-STORM system, from the fuel inside the cask to the external walls of the overpack.

Given that, a set of hypotheses and approximations have been adopted (the first two being scenario-related and the last two being modeling-related):

• The axial heat distribution has been considered uniform;
• The fuel assemblies have been treated as equally powered;
• The bottom/top plug and spacer grids have been not modeled;
• The insertion channels in the down-comer have been ignored (conservative choice).

To model the cask, particular attention has been paid to the creation of two different nodalizations: one to reproduce the thermal response of fuel, the other for the thermal hydraulic behavior of the two gas regions (helium and air). The former is required by the MELCOR COR package and involves radial rings and axial levels, the combination of which results in the creation of cells to correctly represent the “core” region. The latter, instead, involves control volumes (CVs) and flow paths (FLs), addressed through the appropriate parameters to represent the helium flow inside the MPC and the air circulation in the gap between the MPC and the overpack. In particular, friction factors and form losses have been defined carefully, as suggested in [33,34]. In addition, the walls of both the canister and the overpack have been modeled using the MELCOR Heat Structure (HS) package, as in Figure 2. The characteristic lengths for the MPC walls have been imposed equal to their thickness to take into account the main role of the heat exchange from helium to air. Moreover, radiative heat transfer has been managed through the activation of the structure-to-structure model between the honeycomb and the canister walls, and between the canister and the cask. Since the considered surfaces are facing each other and the distance between them is small (17.5 cm in between the honeycomb and the MPC; 4.5 cm between the MPC and the overpack), the view factors have been considered equal to 1.0, meaning that the emissions of a surface are intercepted entirely by the facing one. As heat exchange is a key issue, the definition of radiation parameters and proper material properties has been given utmost importance: in particular, carbon and stainless steel emissivities and concrete conductivity have been user-imposed (as in [22]), whereas other material properties have been taken as in the MELCOR Material Properties (MP) package database.

Finally, the Decay Heat (DCH) package has been activated to impose the decay power (30.0 kW) in the active zone of the fuel assemblies.

3. Uncertainty and Sensitivity Analysis

3.1. General Approach

As said above, PCT is the main thermal indicator of cladding integrity for SNF. Its accurate determination through a realistic but approximate thermal-hydraulic assessment based on UASA methodology (see Section 3.2.1) is the main purpose of this paper. In a first level discussion, the fuel temperature in a dry cask will be dependent on the heat source within the fuel (DCH) and the heat rejection to the sink (environment); the impact of these two factors will be directly explored through a separate sensitivity analysis (SA), as reported in the following. This said, the path and rate at which such heat transfer from the fuel to the environment occurs might also play a role in the final value of the maximum cladding temperature; in this regard, the amount of parameters and complexity of the heat transfer mechanisms involved make it recommendable to use UA to determine the uncertainties associated with the modeling of heat transmission in MELCOR estimates. By considering this, a first-of-a-kind UA is proposed in this work.

3.2. Uncertainty Analysis

3.2.1. Uncertainty Methodology

Several methodologies have been proposed and extensively used in thermal-hydraulics for the evaluation of uncertainties as a support of best-estimate (BE) safety analysis [35,36]. Among these, the GRS methodology [37] has been selected as starting point for the development of a series of MATLAB [38] scripts devoted to the unfolding of the uncertainty analysis carried out in this work.

It is worth highlighting that the chosen methodology is based on the propagation of the input parameters. Given that the only source of uncertainties considered in the analysis is the one related precisely to the selected parameters, the resulting uncertainty value will be smaller than the real one.

The MATLAB scripts follow the logic shown in Figure 3:

• SAMPLING: as the methodology is based on the propagation of input uncertainty parameters, a set of n samples (or variates) is created by means of random sampling (Monte Carlo method) on the basis of user-selected parameters and their distributions. The number of samples is selected according to the minimum number of code calculations required to obtain the probability and confidence levels selected for the UA. According to the Wilks’ formula for one-sided statistical tolerance limits [39], 59 runs are required for a 95%/95% probability and confidence levels (i.e., there is a 95% probability that the maximum value obtained for the selected FOM will not be exceeded);
• INPUT DECK CREATION: the n samples are employed to set up n new MELCOR input decks. Starting from the base case input deck, parameters are identified and replaced with the sampled values;
• BATCH FILE CREATION/EXECUTION: batch files are generated for the execution of the n input decks under the Windows operating system. The batch files contain instructions for the subsequent launch of MELGEN and MELCOR executables, as required by the code;
• STATISTICAL ANALYSIS & CORRELATION COEFFICIENTS: once the results from the calculations were obtained, the user-selected response, hereafter called FOM, was analyzed. For the statistical analysis, the minimum, maximum and mean values as well as standard deviation, are evaluated along the entire calculation time or for the peak value only. As a first approach, the Pearson and Spearman correlation coefficients are used as measures for linear and/or monotonic relationships among each of the input parameters and the FOM [40]. Partial correlations are additionally calculated to take into account the presence of the other input parameters.
• PLOTS: a number of plots are generated, reproducing the results of the uncertainty study, e.g., dispersion plots, Pearson and Spearman visualization plots, probability density functions (PDF) and cumulative distribution function (CDF), FOM vs. iteration index, etc.

3.2.2. Uncertainty Configuration

Given the pioneering application of MELCOR to dry cask modeling in the frame of a BEPU analysis, it was decided to keep the problem as close to the original one as possible and, as a consequence, to not vary the geometry of the system, the masses or the initial and boundary conditions. Input parameters were selected among the ones related to the heat removal capabilities of the canister-overpack system. The choice of parameters has been led by those shown to be more influential [26] and by their accessibility to MELCOR users.

The user-selected parameters and their variation ranges are reported in

Table 4. The input parameters.

Parameter	Base Case	Range of Variation
Stainless Steel Emissivity	0.59	0.1–0.9
Carbon Steel Emissivity	0.66	0.1–0.9
Concrete Conductivity (W/m-K)	1.4	0.7–2.1
FCELR	0.1	0.0–0.5
FCELA	0.1	0.0–0.5
Form Loss Coeff.—Inlet	1.0	0.5–2.0
Form Loss Coeff.—Outlet	1.0	1.0–4.0

. Uniform distributions for these parameters have been chosen in lack of clear data.

Stainless and carbon steel emissivities play a role in the radiative heat exchange between the honeycomb and the stainless steel canister, and between the canister and the carbon steel layer covering the overpack. Not knowing the exact steels that have been employed for the HI-STORM 100S, a broad range of variation has been selected.

Concrete conductivity is the parameter addressing the thermal conduction within the overpack. Given that the overpack is made of plain concrete [25], an upper bound for the conductivity has been set at 2.1 W/m-K. Values higher than this are usually related to reinforced concrete [41,42,43,44].

FCELR and FCELA factors account for the radiative exchange among cells in the “core” region: more precisely, they represent the exchanges occurring radially outward and axially upward, respectively. FCELR and FCELA are “effective” view factors: the default value (0.1) takes into account the fact that, when considering surfaces far from the cell boundary, the line of sight is intercepted by other surfaces, greatly reducing the differential view factors [32]. Values higher than the default one are investigated to evaluate potential effects on the system, also considering that the HI-STORM “core-like” region is smaller than that of a LWR, for which the default value has been derived.

In addition, form losses coefficients for the entrance and exit of the air in the gap between MPC and overpack have been taken into account. When considering sharp entrances or exits, the values for such coefficients can be imposed as 0.5 and 1.0, respectively [45]. Since in the MELCOR model inlet and outlet ducts have been modeled with only one equivalent inlet opening and one equivalent outlet opening, the lower bound for the variation range of the form losses coefficients has been set as in the case of a single inlet and a single outlet. The upper bound, instead, has been set considering four inlets and four outlets, as in the real cask.

It is worth pointing out that some of the parameters related to the heat transfer (e.g., helium properties, view factors between honeycomb and MPC and between MPC and overpack, etc.) have not been considered for the UA. This limitation is mainly linked to the structure of the MELCOR code (and of its input deck), which makes the accessibility (and variation) of some parameters not straightforward.

3.3. Sensitivity Cases

A separate sensitivity analysis has been carried out to assess the influence of initial and boundary conditions on the system, thereby considering scenarios slightly different than the initial one. Two sensitivity cases have been performed in relation to the decay power, which corresponds to the heat source of the system. Five other cases have been run to address the temperature of the air cooling the system, as it acts as the heat sink. The choice of decay power and external air temperature as the main parameters for the sensitivity analysis (SA) is additionally supported by the results obtained in [46], in which both parameters seem to play a role when analyzing the PCT (and the cask external wall) in a CONSTOR RBMK-1500 cask.

For the decay power, the considered values correspond to ±10% of the decay power employed in the base case (30 kW): this choice is supported by [47], where a total uncertainty of 5.6% has been calculated in relation to the decay power of a spent fuel assembly (SFA) of around 50 GWd/tU. As for air, temperatures ranging from 273.15 K to 323.15 K have been taken into account to address the climate in different places in the world.

The sensitivity cases are reported in

Table 5. The sensitivity cases.

-	Cases
Decay Power	27 kW, 33 kW
Air Temperature	273.13 K, 283.15 K, 293.15 K, 308.15 K, 323.15 K

.

4. Results and Discussion

4.1. Base Case Results

Before proceeding with the UASA, a best estimate (BE) case has been calculated. The main BE results obtained are presented here. As PCT is the FOM selected for this work, it is the focus of the results reporting. Additional thermal results might be found in [30].

Figure 4 shows that the PCT (631.35 K) is more than 40 K below the regulatory threshold (673 K). As expected, when DCH distribution is homogeneous across fuel assemblies, the PCT is reached in ring 1 (the most inner radial node of the system meshing). As for axial location, the PCT is in axial level 11, in the upper half of the fuel, which is consistent with the nearly flat profile of DCH axial distribution and the progressive heat-up of helium when moving upward within the MPC.

It should be noted that the time evolution shown in Figure 4 is meaningless, since in reality the cask would reach the steady state from a peak temperature attained during the so called “drying process,” during which the MPC goes through time intervals with no helium supporting the convective heat removal from fuel rods [48].

In an effort to have a sound verification of the obtained results, the MELCOR predictions have been compared to those deriving from a previous CFD modeling of exactly the same system [22]. As reported in

Table 6. The MELCOR-FLUENT comparison—PCT [30].

	PCT (K)	Axial Position (m)
MELCOR	631.35	4.04–4.55 (11th axial level)
FLUENT [22]	628.1	4.22

, the two codes agree on both the PCT value and its location. The small difference between the predicted PCTs is likely within the uncertainty range of the computational tools, and it can be considered acceptable on the account of the very different nature of the employed codes. Besides, the axial location difference comes mostly from the very different meshing schemes used by the two code approximations: lumped parameter codes (MELCOR) versus computational fluid dynamics (FLUENT).

As further confirmation of the MELCOR code capabilities in evaluating the thermal behavior of the dry cask, the temperature radial profile, measured at the height in which the PCT takes place, is considered, as shown in Figure 5. MELCOR provides a discontinuous description of the thermal profile, according to the radial nodalization employed in the “core” zone. Despite this, the predicted thermal behavior can be considered consistent with the CFD one and with the physics behind it. Both codes agree that the PCT is attained at the center of the MPC, which is the region farther away from the cooling effect of the air stream. Fuel temperature decreases when considering locations away from the center, with a small deviation between the codes’ predictions. Deviations become larger in the helium downcomer and in the air channel, which is consistent with the ability of CFD codes to mimic the effect of thermal boundary layers between surfaces and convection fluid. In essence, when comparing the radial thermal profiles predicted by MELCOR and FLUENT, the agreement is substantial. The correspondence between the two codes’ results is remarkable, especially in the fuel zone, but it can be considered acceptable all along the entire profile.

4.2. UASA Outcomes

4.2.1. Uncertainty Analysis

In order to obtain the PCT predictions uncertainty when heat transfer parameters involved in the MELCOR formulation propagate their own uncertainty all the way through the code algorithms, 59 calculations have been successfully run according to the Wilks’ formula for 95%/95% probability and confidence levels. The results obtained through the perturbation of the input uncertainty parameters have then been statistically analyzed with MATLAB.

As can been seen in Figure 6, the dispersion plot shows that the PCT at the steady state does vary less than 15 K, which corresponds approximatively to a discrepancy of 2.5% with respect to the obtained BE value. Considering the order statistics behind the Wilks formula, the maximum value obtained for the PCT corresponds to the 95th percentile with a 95% confidence level for the one-sided tolerance limit. In other words, even if the upper bound of the PCT uncertainty (

Table 7. FOM—statistical analysis.

	Value
PCTmax	637.35 K
PCTmin	623.53 K
PCTmean	630.15 K
std	3.533 K

) is adopted, there is still a nearly 35 K margin to the regulatory limit (673 K).

A supplementary way to present the obtained results is shown in Figure 7, where the probability density function (PDF) and the cumulative distribution function (CDF) have been calculated for the selected FOM. A tentative continuous PDF has been derived considering a normal distribution, as the FOM distribution was not known a priori. The good fit observed indicates that PCT uncertainty distributes rather evenly at both sides of the mean, with no specific PCT shift to higher or lower values and highlights the fact that the PCT BE value is a good representation of the FOM. Namely, despite having adopted default uniform distributions when characterizing the input deck uncertainties, the ranges set were so narrowly that their effect on the PCT is rather moderate and does not drastically change the noticeable safety margin noted for the BE. CDF, instead, provides a graphical visualization of the data obtained from the statistical analysis.

Furthermore, few normality tests [49,50] have been performed to determine if the data set is well-modeled by a normal distribution. Lilliefors, Jarque-Bera, Anderson-Darling and Shapiro-Wilk tests have shown a p-value higher than 0.05, proving that the output distribution can be assimilated to a normal one.

The p-values for the abovementioned normality tests are reported in

Table 8. The p-values for normality tests.

	p-Value
Lilliefors test	0.5
Jarque-Bera test	0.19
Anderson-Darling test	0.42
Shapiro-Wilk test	0.28

.

The influence of each single input parameter on the figure of merit, namely the PCT, has been studied by means of Pearson and Spearman correlation coefficients. More precisely, partial coefficients have been evaluated measuring the intensity of the relationship between each parameter and the FOM, while taking into account the effect of the other parameters.

The results (Figure 8) show no significant relationship, neither linear nor monotonic, between the considered parameters and the PCT. Nonetheless, the results show that in-between metallic surfaces radiation and form losses do play a minor role in the scenario, as does radial radiation within the core region, and not much effort should be devoted to accurately estimating the value of the parameters involved.

4.2.2. Sensitivity Analysis

Concerning the sensitivity analysis, decay power and air temperature have been addressed to assess the behavior of the system following their variation. A total of seven sensitivity cases have been considered: five for the air temperature and two for the decay power. As a result, a linear relationship between the PCT and both air temperature and decay power has been found, as shown in Figure 9 and Figure 10 (in both plots, the black dot indicates the base case). Linear fitting parameters are reported on both plots together with the R² coefficient: in both cases, the R² value is close to 1, indicating a very good fit.

It is interesting to see that a variation of 1 K in the air temperature corresponds to a variation of 1 K (slightly more, to be precise) in the PCT. This could mean a variation of 20 K in the PCT if the cask is located in two sites having a difference of around 20 K between their mean external temperatures. In a more general way, these results indicate that climate has a non-negligible effect on the PCT value. The variation of the PCT with respect to the base case value is reported in Figure 11 for each of the studied cases.

As for the decay power, a 10% variation in the decay power results in a variation of about 25 K in the PCT, corresponding to a percentage change of 4%, as shown in Figure 11 and Figure 12. In this respect, keeping a temperature of 300 K and considering a 10% increase in the decay power (33 kW), the safety margin is reduced to 15 K.

5. Conclusions

The present paper reports a first-of-a-kind uncertainty application for a dry cask system for SNF storage. The MELCOR 2.2 code has been employed as a main computational tool in combination with the MATLAB programming environment for the management of the entire uncertainty and sensitivity analysis.

As for other nuclear systems, USNRC interim staff guidance [7] proposes, for dry casks, acceptance criteria to guarantee fuel rods integrity, with a focus on the maximum value that cladding could reach, namely the PCT. In particular, in [7], a specific threshold temperature of 673 K applicable to all commercial spent fuel burnup levels and cladding materials is given, although for low burnups (less than 45 GWd/tU)a higher short-term temperature limit (843 K) could be accepted.

In this regard, the PCT has been designed as the main FOM. Input parameters for the uncertainty analysis have been selected among the ones related to the heat transfer capabilities of the system, whilst decay power and external air temperature, namely heat source and heat sink, have been chosen as variables for the sensitivity analysis.

The results indicate that the considered parameters have a weak influence on the PCT. On the contrary, an increase in the decay power and in the external air temperature lead to a non-negligible reduction of the safety margin, intended as the difference between the regulatory limit of 673 K and the maximum cladding temperature in the canister.

Notwithstanding the lower computational cost with respect to CFD codes (which are widely used to study dry cask systems), MELCOR seems to have the ability to picture the thermal behavior of the cask. In other words, if a specific safety assessment does not require a full hydro-thermal picture of the dry cask, the study shows that a proper MELCOR model can do an accurate job. This might be of utmost relevance if a full-scope uncertainty and sensitivity analysis of cask estimates is planned to be conducted.

Having said that, further studies could be useful to address more extreme scenarios and to deliver a more comprehensive uncertainty assessment. However, the MELCOR 2.2 code and the developed MATLAB scripts can be considered suitable for the management of the entire uncertainty and sensitivity analysis of a SNF dry cask.