陈仁宗，王 辉，孙晓晖

(1.清华大学 核能与新能源技术研究院，北京 100084；2.中国核电工程有限公司，北京 100840)

1988年，美国核管会(NRC)修订了美国联邦法规10CFR50.46允许使用最佳估算方法进行大破口事故的认证级分析，但必须考虑不确定性，并加以量化计算，以保证分析结果在验收准则之内具有较高的概率[1-4]。最佳估算方法应符合以下3个条件：应根据选定的验收准则，事故分析不引入有意的悲观性；使用最佳估算程序；包含不确定性分析[5]。最佳估算方法利用尽可能详细的模型而非简单模型来保证结果更为接近物理现实，用最佳估算结果的不确定性来量化分析结果与物理现实之间的差距[6]。

目前，最佳估算加不确定性(BEPU)方法主要应用于设计基准事故的安全分析[7-9]。与设计基准事故相比，严重事故现象更加复杂，涉及堆芯裸露升温、燃料包壳氧化升温失效、堆芯熔毁和裂变产物释放等诸多过程。复杂的严重事故现象导致了严重事故更大的不确定性和更为复杂的不确定性分析计算过程。严重事故现象复杂，其不确定性较大，目前尚无官方发布的严重事故现象识别排序(PIRT)表。本文采用的方法是结合设计基准事故PIRT表中参数和与严重事故现象相关的参数进行严重事故输入参数的不确定性分析。

堆芯出口温度(CET)是核电厂安全运行的重要监测参数，在严重事故管理导则(SAMG)中作为堆芯状态的重要表征参数，可作为堆芯损伤评价和简化源项评估的整定值，并可作为严重事故管理导则的入口条件。同时，考虑堆芯裸露后复杂的严重事故现象，文中并未针对包壳破裂后的严重事故进程进行不确定性分析。包壳破裂通常被作为裂变产物气隙释放过程的开始，在堆芯损伤评估和简化源项评估过程中扮演重要的角色。

本文以VVER1000压水堆核电厂为研究对象，选取大破口失水事故(LBLOCA)始发严重事故进行包壳破裂对应堆芯出口温度的BEPU分析。

1 计算模型

本文以VVER1000压水堆核电厂为研究对象，BEPU计算流程如图1所示。第1步，建立VVER1000的最佳估算模型，并对模型进行验证；第2步，选取不确定性分析的输入参数，并对输入参数进行取样，根据取样结果进行计算；第3步，整理计算结果，进行不确定性和敏感性分析。

图1 BEPU计算流程Fig.1 Flowchart of BEPU calculation

本文采用的最佳估算程序为MELCOR。MELCOR是NRC委托桑迪亚国家实验室开发的用于概率风险评估的程序[10]，主要用于严重事故分析，如SAMG中时间窗口的计算、氢气风险和源项评估等，目前也是监管机构主要的安全评审工具。MELCOR的运行需要两步，第1步执行MELGEN进行输入检查，第2步启动MELCOR进行迭代计算。

1.1 VVER1000模型建立与瞬态计算

VVER1000压水堆核电厂是由俄罗斯研发的第3代核电技术[11]。核电厂主回路系统包括4个环路，每个环路有1台主泵和1台直流式蒸汽发生器，共用1台稳压器。反应堆热功率为3 000 MW，堆芯出口处压力为15.7 MPa，入口冷却剂温度为291 ℃。采用MELCOR对VVER1000压水堆核电厂的设计特点建立最佳估算模型，系统节点图如图2所示。建模对象包括反应堆压力容器、稳压器、直流式蒸汽发生器、主泵、主给水系统和主管道等。其中压力容器包括下降段、下封头、下腔室、堆芯区域、旁通段和上腔室的模拟。

在进行瞬态计算前要对模型进行稳态调试，确保模拟值与设计值一致。冷却剂流量、堆芯入口温度、堆芯出口温度、堆芯压力的模拟值与设计值对比如图3所示，模拟值与设计值的相对误差列于表1，由图3和表1可知，模拟值与设计值之间的误差较小，相对误差均在1%以内，说明所建模型可信，可用于进一步的计算。

图2 主回路系统节点图Fig.2 Node diagram of main circuit system

LBLOCA是轻水堆安全设计的基准事故[12]。国际原子能机构把LBLOCA始发严重事故作为计算简化源项的基准事故[13]，国内及美国NRC将NUREG-1465报告中LBLOCA始发严重事故的源项作为核应急的参考源项[14]。本文选取LBLOCA作为始发事故，事故假设条件如下：0 s时刻，反应堆冷管段发生当量直径为30 cm的破口；余热排出系统不可用；安注箱的非能动部分可用，能动部分失效；蒸汽发生器二次侧辅助给水丧失。

事故发生后，大量冷却剂泄漏至安全壳内，导致安全壳内压力升高，同时一回路压力迅速降低，从而触发反应堆紧急停堆及主泵停运。压力降低到一定值时，安注箱非能动部分开始向堆芯注水，随着应急冷却水的不足，堆芯开始裸露升温，然后包壳破裂，开始气隙释放。大破口失水事故序列列于表2。

图3 不同参数的模拟值与设计值对比Fig.3 Comparison of simulation value and design value of different parameters

表1 模拟值与设计值的相对误差Table 1 Relative error between simulation value and design value

表2 事故序列Table 2 Accident sequence

1.2 不确定性分析方法

本文采用拉丁超立方抽样方式对输入参数进行抽样。与简单随机抽样相比，拉丁超立方抽样具有更高的抽样维度，抽取的样品参数更具代表性。拉丁超立方抽样的实现步骤如下：1) 将参数按其范围进行等概率划分，划分为很多子范围；2) 用随机数产生器产生指定数量的[0，1]范围内的随机数；3) 根据均匀分布或正态分布将随机数映射为实际范围内的参数。

本文采用Wilks非参数统计方法[15-16]进行不确定性量化计算，此方法广泛应用于BEPU统计分析中。其基本思想是通过Wilks公式根据选定的置信水平和概率水平获取最小的计算次数，通过对响应参数的有序统计理论来获取统计容忍区间。Wilks公式为：

β=1-γN

(1)

式中：β为置信水平；γ为概率水平；N为计算次数。

基于Wilks公式，对于单侧统计容忍限值(95/95)，样本数决定了限值的选取原则：样本为59组，则许用限值为59组计算结果的最大值；若样本为93组，则许用限值为93组计算结果中的第2大值；若样本为124组，则许用限值为124组计算结果中的第3大值。本文选取的样本数为59。

2 输入参数选取与抽样

2.1 输入参数选取

结合用于设计基准事故的LBLOCA PIRT表[17-18]、与锆包壳破裂相关的严重事故现象和专家判断进行输入参数的确定。综合考虑核电厂设计文件、工程经验判断和相关文献资料中的参考值进行输入参数范围和概率密度分布的判定。不确定性输入参数及其范围、分布列于表3。

表3 不确定性输入参数Table 3 Uncertainty input parameter

2.2 输入参数抽样

输入参数的抽样采用DAKOTA程序[19]实现，选用Mersenne Twister随机数产生器，种子设定值为98 765。对不确定性输入参数分别进行59次拉丁超立方抽样，归一化结果如图4所示。

图4 不确定性输入参数归一化分布Fig.4 Normalized distribution of uncertainty input parameter

由图4可看出，参数的归一化结果与参数分布形式有关，均匀分布参数的归一化结果呈现均匀分布，正态分布参数的归一化结果呈现中间聚集效应。输入参数抽样结果达到了预期的效果。

3 结果分析

3.1 不确定性量化

对59组MELCOR计算结果进行分析整理，得到气隙释放时间和堆芯出口温度，如图5所示。根据Wilks非参数统计方法，对于单侧统计容忍限值(95/95)，样本为59组，则许用限值为59组计算结果的最大值。据此得到的气隙释放时间和堆芯出口温度单侧统计容忍限值分别为1 905.71 s和430.85 ℃。

3.2 敏感性分析

Spearman秩相关系数法是全局敏感性方法，可检验多个输入参数对响应参数的全局影响。采用Spearman秩相关系数法作为输入参数对堆芯出口温度敏感性进行分析，相关系数取值范围为[-1，1]，绝对值的大小表示敏感性的强弱；数值的正负号表示正负相关性，正号表示正相关性，负号表示负相关性，零表示不相关，其表达式如式(2)所示。

图5 气隙释放时间和堆芯出口温度Fig.5 Gas release time and CET

ρ=

(2)

式中：RXi为Xi在X中的大小排序；RYi为Yi在Y中的大小排序；ρ为秩相关系数；n为样本数。

输入参数与堆芯出口温度间的Spearman秩相关系数如图6所示，计算结果表明，气隙释放对应的堆芯出口温度对衰变热系数和包壳厚度较为敏感。

图6 输入参数与堆芯出口温度的 Spearman秩相关系数Fig.6 Spearman rank relational coefficient of input parameter and CET

4 结论

以VVER1000压水堆核电厂LBLOCA始发严重事故为研究对象，MELCOR为研究工具，开展了输入参数对堆芯出口温度的不确定性和敏感性研究，得到的结论如下：1) 气隙释放时间和堆芯出口温度单侧统计容忍限值(95/95)分别为1 905.71 s和430.85 ℃；2) 气隙释放对应的堆芯出口温度对衰变热系数和包壳厚度较为敏感。