Abstract
Uncertainty of a severe accident code output needs to be handled reliably considering its use in safety regulation of nuclear industry. In particular, severe accident codes are utilized for probabilistic safety assessment (PSA), where the uncertainty of severe accident progress should be considered carefully due to its influence on human reliability analysis. Therefore, in this study, the uncertainty analysis of severe accident progress was performed using MELCOR code, and a total of 200 data sets of in-vessel uncertainty parameters were generated by Latin hypercube sampling method. The rank regression analysis was also performed to investigate the effect of uncertainty parameters on the severe accident progress. Sensitivity coefficients (SCs) in MELCOR such as molten clad drainage rate and zircaloy melt breakout temperature showed significant influence on relocation time and dryout time of lower plenum. However, the influence of uncertainty parameter diminished as the accident progressed.
1. Introduction
Since Fukushima accident in 2011, numerous studies about the severe accident (SA) phenomena and their progress have been investigated experimentally and numerically. However, due to the large scale of nuclear power plants, not many experimental studies were conducted for the accident sequence. Therefore the researches about the SA progress were conducted mostly by the numerical analysis [1–5]. As far as the numerical analysis of SA progress is concerned, lumped parameter SA codes such as MELCOR [6], MAAP [7], and SCDAP/RELAP [8] have been utilized in general. However, it is widely accepted that the SA code results may generate large uncertainty because rather simplified models were adopted in the codes due to the incomplete knowledge about SA phenomena. Because the code results can be reutilized for calculating source terms in the frame of PSA level 3, it is very important to investigate the uncertainty range and to identify relative importance of input parameters [9, 10].
Among the code outputs, identifying timing of certain SA sequence is important in terms of human reliability analysis (HRA). As a practical example, the success probability of human related actions can be defined as a function of available time in HRA [11]. Thus, the timing of SA progress is important to define the available time for certain human related actions. Therefore, the uncertainty of the SA progress in terms of the major sequence time needs to be analysed for better understanding on the nuclear safety. For this reason, the major objective of this study is to perform the uncertainty analysis about occurrence time of SA progression phenomena. A representative SA scenario was selected as short term station black out (STSBO). STSBO is the scenario that all power, even direct current (DC) batteries, is unavailable. Because the safety feature using batteries is also unavailable, the accident progress can be faster than other accident scenarios.
2. Numerical Methodologies
2.1. Description of MELCOR Code [6]
MELCOR code is an integrated SA code for various kinds of nuclear power reactors. Since 1982, MELCOR code has been developed by Sandia National Laboratories (SNL) and used for plant risk assessment and source term analysis. To simulate various SA phenomena, a number of modular packages were coupled in a unified frame. By the interaction of each modular package, MELCOR code can simulate thermal-hydraulic response of the plant, heat-up, oxidation, degradation and relocation of the core, and fission product release and transport. Because of the lack of knowledge about the SA phenomena and the interest of quick code execution time, parametric equations were used to model complicated physical processes in the past. However most recent MELCOR models are mechanistic thank to the improved calculation speed of modern computers. In many of the mechanistic models, adjustable optional parameters can be also available for uncertainty analyses and sensitivity studies. For that reason, MELCOR version 2.1.6342 was used.
2.2. MELCOR Input Model of OPR1000
As a reference plant, Optimized Power Reactor 1000 MWe (OPR1000), which consists of majority of Korean operating NPP, was selected. The nodalization of OPR1000 input model is shown in Figure 1. The input model includes reactor coolant system (RCS) composed of 25 control volumes (CVs) and simple secondary side composed of two steam generators (SGs) and a turbine. RCS consists of two hot legs (CV310 and 410), four cold legs (CV380, 390, 480 and 490), pressurizer (CV500), and reactor pressure vessel (RPV) with 6 detailed volumes (CV130, 150, 180, 170, 190, and 260). At the top of a pressurizer, pressurizer safety relief valves (PRSVs) and safety depressurization system (SDS) were modeled to simulate release of RCS coolant in case of depressurization. The containment (CNMT) consists of 20 CVs and the total free volume is 77,440 m3. Figure 2 shows the nodalization of the core. Total 7 ~ 8 radial rings and 14 axial nodes were modeled for the detailed simulation of the reactor core. 1st ~ 3rd axial core nodes were coupled with CV150 (LP), and 4th ~ 14th axial core nodes were coupled with CV170 (core region). In the core region, only 4th ~ 13th axial core nodes contained the nuclear fuel. The supporting structures were located at the 1st, 2nd, 3rd, 6th, 9th, and 12th axial node, and support plate was modeled in 3rd axial node. Total initial mass of uranium fuel, zirconium cladding, and stainless steel in the core were 85,692.0 kg, 23,910.58 kg, and 20,024.4 kg, respectively.
2.3. Description of Uncertainty Parameters and Distribution
In SOARCA project [9], the uncertainty analysis of source term was performed using MELCOR code. Based on PIRT process by SNL and nuclear regulatory commission (NRC), total 24 epistemic MELCOR parameters were selected for uncertain input parameters of Surry power plant. To study the effect of molten corium behavior on the severe accident progress, 5 of uncertain parameters in COR package were selected in this study. The detailed distribution of uncertain parameters referred to the uncertainty analysis of Surry plant in SOARCA project. Table 1 shows the selected uncertain parameters and distributions. The lower/upper bound of zircaloy melt breakout temperature (SC1131(2)) was defined by melting point of zircaloy (2100 K) and the value of hydrogen uncertainty study [12] based on Phebus experiments [13]. The distribution was updated from the distribution of Peach bottom plants with the S/Q simulation results [14]. Because the technical base for the boundary values of molten clad drainage rate (SC1141(2)) and radial solid/molten debris relocation time constant (SC1020) were insufficient, the boundaries of these uncertainty parameters were selected based on order of magnitude. The default value of molten clad drainage rate was determined with CORA-13 experiment. Therefore, log triangular distribution, that is, the mode value is this default value, was used for the distribution of molten clad drainage rate. The uniform distribution was used for radial solid/molten debris relocation time constant due to the lack of experimental data. The distribution of Effective temperature at which the eutectic formed from UO2 and ZrO2 melts was created based on the mean and standard deviation of 6 test results from VERCORS experiment [15].
2.4. Description of Sampling and Quantification of Results
Symbolic nuclear analysis package (SNAP) is a graphical user interface for simplifying the task of creating input files for the analytic codes and helping to visualize code results [16]. Input files of MELCOR code can also be created by SNAP. For uncertainty and sensitivity analysis, SNAP provides the additional plugin called “DAKOTA”. Using DAKOTA plugin, users can create the random samples of each sensitivity coefficient by LHS method and combine them into an input data set. The correlation between uncertain parameters was not considered in this study. Variance inflation factor and correlation coefficient between sampled input parameters were checked to prevent significant correlation between input parameters by random sampling. In this study, total 200 cases were generated and simulated until 259,200 s (72 hrs). For the quantification of the effect of the selected uncertain parameters on the accident progress, the time of initial melt relocation, LP dryout, RPV failure, and CNMT leak were investigated as major output parameters. In addition, the statistics of the major output parameter were calculated. To investigate the individual effect of uncertain input parameters, linear rank regression analysis was also performed.
3. Results and Discussion
3.1. Description of Base Case
To investigate the uncertainty of SA progress, the aforementioned STSBO was selected as a base case. All of safety features were assumed to fail by losing AC power, and steam driven auxiliary feed water (AFW) pump was also modeled to fail. Table 2 shows the SA sequence of the base case. The STSBO was initiated by losing AC power and subsequently the reactor was tripped. With the reactor tripped, reactor coolant pumps (RCPs) and main feed water (MFW) pumps stopped. Following the reactor trip, the decay heat from the core was removed by the natural circulation of primary coolant during early phase of the accident. Because the AFW pumps were unavailable by the assumption, all SGs were exhausted at 3,762 s by the heat from the RCS. The RCS pressure increased with the loss of secondary heat removal, and the PSRV was opened at 4,881 s as shown in Figure 3. The reactor coolant was released to the CNMT through the PSRV, and thus it caused the core uncovery at 5,359 s. The exposed core was heated up due to the lack of coolant, and additional oxidation heat generated between high temperature steam and high cladding temperature (>1,100 K) accelerated the core heat-up from 8,991 s. The coolant in the core region was dried out at 9,990 s and dryout of the LP region occurred at 12,866 s due to the relocated core debris as shown in Figure 4. After the first LP dryout, thru-wall yielding mode of RPV failure was initiated by the heat of relocated core debris.
When the RPV failure was initiated, the RCS pressure reached to 16.25 MPa (Figure 3). Due to the high pressure, the core debris was released into the CNMT via high pressure melt ejection (HPME) phenomena. Figure 5 shows the ejected debris mass from the RPV, debris mass in the cavity, and mass of debris ejected to CNMT atmosphere by the HPME. At first ejection with RPV failure (14,462 s), 6,341 kg of the 8,819 kg of debris ejected from the RPV was ejected to CNMT atmosphere. Because the RCS pressure decreased after the first ejection, the ejected debris mass following the first ejection only contributes to the increase of debris mass retained in the cavity. After about 19,200 s, 119,870 kg of core debris was ejected to the cavity due to the collapse of remained fuel structure. Figure 6 shows the CNMT pressure during the STSBO accident. Before the RPV failure, the CNMT pressure was repeated to increase and decrease by the opening and closure of the PSRV. After the RPV failure, 0.26 MPa of pressure peak was observed due to the direct containment heating (DCH) by the HPME. Because the debris within the cavity released the heat and noncondensable gas by molten corium concrete interaction (MCCI), the CNMT pressure increased continuously and reached to design pressure (0.44 MPa) at 60,100 s. The CNMT pressure at the end of calculation (259,200 s) was calculated as 1.17 MPa.
3.2. Statistics of Major Accident Progress
Among a total of 200 sample cases, only 162 cases were calculated until 259,000 s (72 hrs). Because of COR or CAV package errors, 37 cases were stopped as soon as the RPV failure was reported. A case was calculated until the initial melt relocation. However, it was terminated before the RPV failure. For that reason, the statistics of accident progress such as mean, coefficient of variance (C.V.), median, 5th percentile (P5), and 95th percentile (P95) were calculated through the available cases. Table 3 shows the statistical values of accident progress. The selected uncertainty parameters affected the accident progress after cladding melt. Therefore, the statistical values of initial melt relocation, first LP dryout, RPV failure, and CNMT overpressurization time were investigated. The mean values of initial melt relocation, first LP dryout, RPV failure, and CNMT overpressurization were calculated as 10,445 s, 13,421 s, 15,850 s, and 60,100 s, respectively. All mean values of accident progress were larger than those of the base case. The values of C.V. show that the variance of accident progress increased as the accident progressed. The median values of initial melt relocation, first LP dryout, RPV failure, and CNMT overpressurization were estimated as 10,491 s, 13,466 s, 15,731 s, and 58,932 s, respectively. Median values of initial melt relocation and first LP dryout time were larger than mean values of those. It means that the distributions were biased toward the high values. On the other hand, RPV failure and CNMT overpressurization time were biased toward the low values. Figure 7 shows these biased distributions of accident progress. Counts on Y-axis corresponded to the number of samples that was included in specific range about the time of accident progression phenomena on X-axis. The 5th/95th values of initial melt relocation, first LP dryout, RPV failure, and CNMT overpressurization were 10,244 s/10,601 s, 12,831 s/13,811 s, 14,911 s/17,252 s, and 55,882 s/69,272 s, respectively. Using bias-corrected and accelerated bootstrap method, 95% confidence interval of 5th/95th values was also calculated. The confidence intervals for 5th/95th values of initial melt relocation, first LP dryout, RPV failure, and CNMT overpressurization were 10,236 ~ 10,246 s/10,586 ~ 10,651 s, 12,676 ~ 12,931 s/13,781 ~ 13,841 s, 14,826 ~ 14,973 s/17,192 ~ 17,304 s, and 55,621 ~ 56,187 s/68,798 ~ 69.639 s, respectively.
(a) Relocation time
(b) First LP dryout time
(c) RPV failure time
(d) CNMT overpressurization time
3.3. Regression Analysis of Uncertainty Parameters on Accident Progress
Although the ordinary linear regression can show the potential relations between input and output values, capturing any complex relationship is challenging in reality [9]. Because the accident progress is the cumulative result of various phenomena, the relationship between in-vessel uncertain parameters and accident progress may be nonlinear or complex. Thus, in this study, rank regression technique was used to investigate the potential relations between in-vessel uncertainty and accident progress. The rank regression is a kind of linear regression, which uses rank values rather than original values. Using the rank values, the rank regression can capture the potential relations including nonlinear influence [17, 18]. For the rank regression results, standardized rank regression coefficient (SRRC) and partial rank correlation coefficient (PRCC) were summarized in Figure 8.
(a) Relocation time (Adj. R2 = 0.612)
(b) First LP dryout time (Adj. R2 = 0.467)
(c) RPV failure time (Adj. R2 = 0.125)
(d) CNMT overpressurization time (Adj. R2 = 0.186)
The coefficient of determination values shows how the regression model explains the relationship between input and output parameters. In this study, adjustable R square (Adj. R2) was used for the coefficient of determination. The Adj. R2 values of the regression model about initial melt relocations show that the in-vessel uncertainty parameters explain the potential relationships about initial melt relocation in a limited way. However, as the accident progressed, the influence of in-vessel uncertainty parameters on the regression model decreased. In cases of RPV failure and CNMT overpressurization, Adj. R2 values were significantly low (Adj. R2 < 0.3). Thus the regression model could not explain the relationship between input and output. These results might be caused by the complex behavior of molten corium relocation. In this case, although the uncertainty parameters affected the tendency of molten corium behavior, the effect on the results can be insignificant due to the cumulative random variation.
In case of initial melt relocation, SRRC and PRCC show strong negative correlation (r < -0.7) between SC1141(2) and the timing of initial melt relocation. Figure 9 shows the scatter plot of relocation time versus SC1141(2). The distribution was divided into two categories by 10,400 s of the relocation time. The strong negative correlation resulted from the upper relocation time category. More cases in the upper relocation time category were observed in lower value of SC1141(2) due to the blockage effect. During the candling process, the molten core materials can be refrozen by the candling heat transfer and form the blockage. Because the lower value of SC1141(2) indicates the slower downward flow of molten core materials, the possibility of blockage formation increased with the lower value of SC1141(2). Therefore, the increment of blockage could delay the relocation time.
SRRC of first LP dryout shows that SC1141(2) and SC1131(2) were the most significant parameters affecting the first LP dryout. SRRC and PRCC indicated moderate negative relationship (-0.3 > r > -0.7) between these two uncertainty parameters and first LP dryout time. Figure 10 shows the scatter plot of first LP dryout time versus SC1131(2). Using the color scatter, the values of SC1141(2) were shown in the scatter plot. The linear fitting line shows the negative relationship between SC1131(2) and first LP dryout time. In case of SC1141(2), more blue color scatters were observed below the fitting line. However, this trend was not clear due to the number of samples below the fitting line. The hydrogen production mass and relocated mass before first LP dryout can describe the effect of SC1131(2) and SC1141(2) on first LP dryout time. Figure 11 shows the SRRC and PRCC values of in-vessel hydrogen production and relocated mass before first LP dryout. The high Adj. R2 value (Adj. R2 > 0.7) shows that the rank regression model can properly explain the relationship between uncertainty parameters and hydrogen production and relocated mass. PRCC indicated strong negative correlation between SC1141(2) and in-vessel hydrogen production and strong positive correlation between SC1141(2) and relocated mass before first LP dryout. While the lower hydrogen production reduced the total heat accumulated in the core, the higher relocated mass transferred more heat to water in LP. Therefore, the effect of SC1141(2) on first LP dryout can be reduced by the countering effect. In case of SC1131(2), only moderate positive correlation between SC1131(2) and in-vessel hydrogen production was observed. For that reason, the effect of SC1131(2) on in-vessel hydrogen production directly influenced the first LP dryout time.
(a) In-vessel hydrogen production (Adj. R2 = 0.704)
(b) Relocated mass before first LP dryout (Adj. R2 = 0.756)
4. Conclusions
To investigate the uncertainty of accident progress induced by in-vessel uncertainty parameters, uncertainty quantification was performed with regard to the timing of initial melt relocation, first LP dryout, RPV failure, and CNMT overpressurization. A total of 200 data sets of sensitivity coefficients in COR package were generated by SNAP-DAKOTA plugin using the LHS method. 5th/95th range of initial melt relocation, first LP dryout, RPV failure, and CNMT overpressurization were 10,244 s/ 10,601 s, 12,831 s/13,811 s, 14,911 s/17,253 s, and 55,882 s/69,272 s, respectively. The rank regression analysis was also performed to investigate the effect of in-vessel uncertainty parameters. The regression results show that zircaloy melt breakout temperature (SC1131(2)) and molten clad drainage rate (SC1141(2)) exhibited significant influence on initial melt relocation time and first LP dryout time, respectively. Although only in-vessel uncertainty parameters among input parameters were modified, the effect on RPV failure and CNMT overpressurization time could not be described by the rank regression model. To clearly understand the effect of in-vessel uncertainty on accident progress, various regression techniques and the stochastic studies about the core behavior should be performed for the future work.
Nomenclature
| Adj. R2: | Adjustable R square |
| AFW: | Auxiliary feed water |
| C.V.: | Coefficient of variance |
| CNMT: | Containment |
| CV: | Control volume |
| DCH: | Direct containment heating |
| HPME: | High pressure melt ejection |
| HRA: | Human reliability analysis |
| LHS: | Latin hypercube sampling |
| LP: | Lower plenum |
| MCCI: | Molten corium concrete interaction |
| MFW: | Main feed water |
| P5: | 5th percentile |
| P95: | 95th percentile |
| PDF: | Probability distribution function |
| PRCC: | Partial rank correlation coefficient |
| PSA: | Probabilistic safety assessment |
| PSRV: | Pressurizer safety relief valve |
| r: | Value of coefficient |
| RCS: | Reactor coolant system |
| RCP: | Reactor coolant pump |
| RPV: | Reactor pressure vessel |
| SA: | Severe accident |
| SC: | Sensitivity coefficient |
| SDS: | Safety depressurization system |
| SG: | Steam generator |
| SANP: | Symbolic nuclear analysis package |
| SNL: | Sandia national laboratories |
| SRRC: | Standardized rank regression coefficient |
| STSBO: | Short term station black out. |
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was supported by the National Research Foundation of Korea (NRF) funded by the Ministry of Science & ICT (MSIT, Korea) [Grant no. NRF-2017M2A8A4018213] and the Nuclear Safety Research Program through the Korea Foundation of Nuclear Safety (KoFONS) using the financial resource granted by the Nuclear Safety and Security Commission (NSSC) of the Republic of Korea (no. 1805001).