1 引言
压水堆小破口失水事故(SBLOCA)属于冷却剂装量减少事故，该类事故因其导致堆芯熔化的风险较高而备受关注。安全评估过程需要对SBLOCA进行全面而深入的分析[1]，为应对各类不确定性，早期常采用保守性安全评价方法，该方法常采用刻意的悲观和简化模型、保守的边界，但这种评价方法计算得到的结果往往过于保守，鉴于该方法的局限性，1988年，美国核管会NRC对法规10 CFR 50.46进行了修订提出可以采用最佳估算方法，但最佳估算的计算结果必须加上不确定性量化评估，即BEPU(best-estimate plus uncertainty)方法，确保计算结果的可信度。
近年来，BEPU方法在世界范围内得到广泛使用，逐步成为核反应堆开展安全评价的一种趋势。2005年，国际原子能机构IAEA同意将BEPU方法用于确定取证所需的安全裕度[2]。最佳估算的不确定性分析方法包括：① 输入不确定性的传播，包括源系数以及输入参数的不确定性传播；② 输出不确定性的外推。本文中利用基于机理模型的最佳估算程序RELAP5/SCDAPSIM对小破口失水事故进行瞬态过程模拟计算，对源系数所带来的不确定性传播进行量化评估，并对计算结果进行敏感性分析，为小型压水堆的安全分析提供理论支撑依据。
2 不确定性量化评估
最佳估算模型安全评估方法是建立在一整套物理模型和相互关联的基础上，在对失水事故进行安全分析时，完整性不仅需要包括流体力学、传热、反应堆动力学、仪表和控制(I&C)等一大套模型，还包括燃料棒性能和安全壳响应模型[3]。
2.1 事故的选取
小型压水堆发生小破口失水事故时，一回路系统压力下降，冷却剂由过冷状态转变为饱和喷放阶段；冷却剂随着破口不断流失，稳压器液位下降、达到高压安注整定值则高压安注自动投入；若安注投入的冷却剂流量比破口处流失的冷却剂流量小，压力继续下降，当压力下降到停堆保护设定值时，反应堆停堆；压力进一步下降，达到低压安注整定值，低压安注自动投入，该过程力图维持堆芯淹没，实现堆芯淹没即可保证堆芯的安全。小破口失水事故界定反应堆是否发生堆芯损毁的依据是燃料包壳是否发生破损，即第一道安全屏障的破损，一般以燃料包壳峰值温度PCT作为评价标准[4]，输出结果重点关注燃料包壳峰值温度。安全分析报告结果表明双环路压水堆破口发生在稳压器侧后果更为严重。
因此，本文采用最佳估算程序RELAP5/SCDAPSIM进行建模，对安全注射系统采用单一故障准则，仅有一台高压安注泵和一台低压安注泵投入，对稳压器侧回路冷管段的 10 mm小破口瞬态过程进行模拟计算。
2.2 现象识别和排序表(PIRT)
基于最佳估算程序RELAP5/SCDAPSIM的事故分析，程序本身采用的机理模型以及输入参数的填写都会带来大量的不确定性传播，本文中将重点针对机理模型开展不确定性分析。
最佳估算程序中与传热相关的不确定性源系数有：① 相间界面传热系数；② 壁面传热系数；③ 临界热通量；④ 来自气隙导热模型的气隙热导率；⑤ 粘性；⑥ 冷却剂导热系数；⑦ 表面张力。根据小破口失水事故的主要现象以及工程判断，对程序中自带的源系数进行一个初步筛选及分级，如表1所示。针对表中所列7种源系数开展不确定性分析，验证本文中的推测。
表1 不确定性源系数分级表
Table 1 The PIRT of uncertainty source code coefficients
影响程度源系数高↓低相间界面传热系数壁面传热系数冷却剂导热系数粘性表面张力临界热通量来自气隙导热模型的气隙热导率
2.3 源系数概率分布及范围
不确定性分析首先需要了解不确定性源系数的概率分布及其特征参数。目前针对小型压水堆的SBLOCA不确定性分析应用较少，缺少相关数据及信息，对于不确定性源系数概率分布及范围参考已发表的相关研究文献[5-6]。本文中选择的不确定性源系数的概率分布及范围见表2所示。
2.4 不确定性分析方法
不确定性分析常采用蒙特卡罗方法，该方法需要进行大量的采样。样本采样主要有2种方法：① 参数抽样；② 非参数抽样。本文中不确定性分析计算所采用的方法是基于Wilks非参数统计的简单随机抽样，样本容量仅由Wilks公式计算得到，与输入参数的个数无关，不需要了解总体分布函数，适用面广。
表2 不确定性源系数的概率分布及范围
Table 2 The probability distribution and range of uncertainty source code coefficients
编号源系数类型描述分布类型注释范围01来自气隙导热模型的气隙热导率气隙导热均匀分布气隙导热模型(0.8,1.2)0203040506070809101112131415相间界面传热系数泡状流弹状流环形雾状流临界热流前雾状流雾状流临界热流后雾状流水平分层流均匀分布液体侧传热(0.8,1.2)均匀分布蒸汽侧传热(0.8,1.2)均匀分布液体侧传热(0.8,1.2)均匀分布蒸汽侧传热(0.8,1.2)均匀分布液体侧传热(0.8,1.2)均匀分布蒸汽侧传热(0.8,1.2)均匀分布液体侧传热(0.8,1.2)均匀分布蒸汽侧传热(0.8,1.2)均匀分布液体侧传热(0.8,1.2)均匀分布蒸汽侧传热(0.8,1.2)均匀分布液体侧传热(0.8,1.2)均匀分布蒸汽侧传热(0.8,1.2)均匀分布液体侧传热(0.8,1.2)均匀分布蒸汽侧传热(0.8,1.2)1617181920212223壁面传热系数单相液体单相气体泡核沸腾过渡沸腾膜沸腾均匀分布(Mode2)(0.8,1.2)均匀分布(Mode9)(0.8,1.2)均匀分布过冷态(Mode3)(0.8,1.2)梯形分布饱和态(Mode4)[0.76,0.86,1.19,1.43]梯形分布过冷态(Mode5)[0.76,0.86,1.19,1.43]梯形分布饱和态(Mode6)[0.76,0.86,1.19,1.43]梯形分布过冷态(Mode7)[0.76,0.86,1.19,1.43]梯形分布饱和态(Mode8)[0.76,0.86,1.19,1.43]24临界热通量CHF查表均匀分布乘数到表(0.8,1.2)25粘性液态轻水粘性正态分布-(1.0,0.02)26冷却剂导热系数液态轻水导热正态分布-(1.0,0.02)27表面张力轻水正态分布-(1.0,0.02)
Wilks方法是一种有序的容忍限值方法，其数学思想可描述如下：对任意的变量X，为其建立一个容忍置信区间(L,U)，使得X的所有取值至少有α份额落在此置信区间的置信度为β[7]。Wilks经典公式如下：
单侧容忍区间：
β=1-αN
(1)
双侧容忍区间：
β=1-αN-N(1-α)αN-1
(2)
式中： β为置信度；α为概率水平；N为样本容量。
美国核管会NRC认为95%置信度，95%的概率水平足以满足安全准则要求，即“95/95准则”。见表3所示，在单侧容忍区间下，当置信度β=95%，概率水平α=95%，程序最少成功运行59次即可满足安全分析要求；同样情况下，在双侧容忍区间，程序最少成功运行93次。程序使用Wilks公式计算所需的不确定性运行次数，代码给定：① 容许限值百分比；② 置信水平；③ Wilks公式应用顺序。程序运行的次数可由上述信息经计算得到，也可由用户自定义。百分位数和置信水平是要获得的不确定性带的特征参数：百分位数是在具有一定置信水平的不确定性界限之下包含的总体数量。用户应指定百分位数和置信水平，对于相应的百分位数和相同的置信水平，还将获得上限和下限。经过3个阶段后，在过程结束时获得的这2个界限被称为单边容许水平，必须分开考虑。
表3 单侧和双侧统计容许极限的最小计算次数N
Table 3 Minimum number of calculations of allowable limits for unilateral and bilateral statistics
单侧容忍区间0.900.950.99双侧容忍区间0.900.950.990.90224523038773880.95295929946934730.99449045964130662
Wilks公式的形式是不完全β函数。代码通过迭代过程计算所需的运行次数，该迭代过程测试“设置”阶段输入卡中提供的置信水平是否大于或等于代码使用不完整β函数计算的置信水平。Wilks公式应用的顺序增加了所需代码运行的次数，但也将获得更准确的不确定性界限估计。
3 结果与讨论
源系数及输入参数的不确定性量化评估可以根据统计参数的分布特点确定每个参数的分布函数类型及范围，相关的分布类型主要有：均匀分布、正态分布、对数—正态分布、梯形分布，常用的分布类型为均匀分布和正态分布。
3.1 正态分布与均匀分布的差异
表2中源系数涉及了2种常用的分布函数，即正态分布和均匀分布，具体采用何种分布函数，通常根据经验选取。本文中为了验证正态分布与均匀分布的区别，以相间界面传热系数为例，采用单一变量法，正态分布范围取(1.0,0.02)，均匀分布范围取(0.8,1.2)，分析相间界面传热系数的不确定性对稳压器压力的影响。从图1、图2的结果中可以发现，正态分布上下限最大跨度为0.303 MPa，均匀分布上下限最大跨度为1.69 MPa，均匀分布相对正态分布的结果不确定性更大，不确定性分析的结果更加保守。
图1 相间界面传热系数为正态分布对稳压器压力的影响曲线
Fig.1 The influence of normal distribution of interphase heat transfer coefficient on pressure of voltage regulator
图2 相间界面传热系数为均匀分布对稳压器压力的影响曲线
Fig.2 The influence of uniform distribution of interphase heat transfer coefficient on pressure of voltage regulator
3.2 不确定性计算方法分析
程序模拟计算包括基本案例运行和设置的由Wilks公式计算所得的59次模拟运行，运行中出现1次计算失败的情况，额外增加1组模拟运行，保证程序最少成功运行59次，满足“PCT95/95准则”[8]。在“模拟”阶段生成的输出值根据其在每个时间步长从低阶(1阶)到高阶(N阶)的等级进行排序，根据用户定义的第一个顺序，代码将使用第一个和最后一个N顺序。经过后处理阶段的计算处理，代码将生成包含以下时间历史的图表：(a)不确定性上限和下限；(b)基本案例时间趋势；(c)每个时间步长上下限之间的跨度。
计算结束，选择小破口失水事故几个关键的输出变量对不确定性的计算结果影响进行分析讨论。图3所示为采用程序自带不确定计算包计算的燃料包壳峰值温度随时间的变化。本文中经过59组基于Wilks公式非参数抽样统计计算得到“PCT95/95”为1.23(归一化处理)，小于1 477 K，满足热工准则要求，可认为此事故条件下反应堆是安全的。不确定带宽最大为82.64 K，发生在停堆瞬变过程，平均值为5.75 K，总体看源系数的不确定性影响相对较小。
图3 燃料包壳峰值温度不确定带曲线
Fig.3 Peak fuel cladding temperature uncertain bandwidth
3.3 敏感性分析
本文中仿真模拟每一种源系数所带来的不确定性传播对PCT的影响(见图4～图10)，通过敏感性分析可以发现，对PCT的影响程度从高到低可依次将不确定性源系数重新排序为：① 壁面传热系数(上下限跨度最大为52.48 K，平均值为2.16 K)；② 冷却剂导热系数(上下限跨度最大为46.39 K，平均值为0.83 K)；③ 相间界面传热系数(上下限跨度最大为7.51 K，平均值为3.23 K)；④ 粘性(上下限跨度最大为0.97 K，平均值为0.38 K)；⑤ 表面张力(上下限跨度最大为0.65 K，平均值为0.18 K)；⑥ 临界热流模型(上下限跨度最大为0.41 K，平均值为0.11 K)；⑦ 气隙导热模型(上下限跨度最大为0 K，平均值为0 K)。
图4 壁面传热系数对PCT的影响曲线
Fig.4 The influence of wall heat transfer coefficient on PCT
图5 冷却剂导热系数对PCT的影响曲线
Fig.5 The influence of thermal conductivity of coolant on PCT
图6 相间界面传热系数对PCT的影响曲线
Fig.6 The influence of interphase heat transfer coefficient on PCT
图7 粘性对PCT的影响曲线
Fig.7 The influence of viscosity on PCT
图8 表面张力对PCT的影响曲线
Fig.8 The influence of surface tension on PCT
图9 临界热流对PCT的影响曲线
Fig.9 The influence of critical heat flux on PCT
图10 气隙导热对PCT的影响曲线
Fig.10 The influence of air gap heat conduction on PCT
综上，对PCT有重大影响的关键不确定性源系数从高到低依次是：壁面传热系数、冷却剂导热系数和相间界面传热系数。可见，针对本文中选取的7种源系数，进行小破口失水事故的不确定性计算时主要考虑上述3种源系数，而其他4种源系数对小破口失水事故的计算没有带来不确定性传播，可以不考虑。
4 结论
以小型压水堆为研究对象，基于最佳估算程序RELAP5/SCDAPSIM对小破口失水事故开展不确定性方法研究以及敏感性分析，主要结论如下：
1) 研究了正态分布与均匀分布的区别，结果表明均匀分布相对正态分布的结果不确定性更大，不确定性分析的结果更加保守；
2) 不确定性计算的结果表明“PCT95/95”为1.23(归一化处理)，远小于热工准则要求的1 477 K，表明在给定的假设条件下反应堆是安全的；
3) 小破口失水事故的安全特性分析，应重点关注对PCT有重大影响的关键不确定性源系数，从高到低依次为：壁面传热系数、冷却剂导热系数和相间界面传热系数。
4) 计算涉及的7种源系数不确定性传播影响不显著，表明其适用于小型压水堆的安全分析计算。
[1] 朱继洲.核反应堆安全分析[M].西安:西安交通大学出版社,2004:1-6.
Zhu J Z.Nuclear reactor safety analysis[M].Xi’an:Xi’an Jiaotong University Publisher,2004:1-6.
[2] 陈炼,房芳芳,邓程程,等.核电站最佳估算安全分析中的不确定度评估方法分析[J].原子能科学技术,2015,49(07):1237-1242.
Chen L,Fang F F,Deng C C,et al.Analysis of uncertainty evaluation method in best e-stimate safety analysis of nuclear power plant[J].Atomic Energy Science and Technology,2015,49(07):1237-1242.
[3] Martin R P,Alessandro P.Progress in international best estimate plus uncertainnty analysis methodologies[J].Nuclear Engineering and Design,2021,374:111033.
[4] RIICHIRO OKAWA,MASAHIRO FURUYA.Large-break LOCA analysis with modified boiling heat-transfer model in TRACE code[J].Nuclear Engineering and Design,2019,346:97-111.
[5] 邓程程,常华健,陈炼.小破口失水事故最佳估算的不确定性和敏感性分析[J].原子能科学计算,2016,50(07):1224-1231.
Deng C C,Chang H J,Chen L.Uncertainty and sensitivity analysis of SBLOCA best-estimate calculation[J].Atomic energy science and technology,2016,50(07):1224-1231.
[6] Yang Y, Yang J, Deng C C,et al.Simulation and uncertainty analysis of main steam line break accident on an integral test facilty[J].Annals of Nuclear Energy,2020,144:107565.
[7] 陈炼,胡啸,邓程程,等.最佳估算方法中不确定度评估关键问题分析[J].原子能科学技术,2016,50(05):851-858.
Chen L,Hu X,Deng C C,et al.Analysis on key issue of uncertainty evaluation in best estimate method[J]. Atomic Energy Science and Technology,2016,50(05):851-858.
[8] 张中伟.保守分析与最佳估算相结合之LOCA认证分析方法论[D].上海:上海交通大学,2011.
Zhang Z W.Deterministic-Realistic Hybrid Methodology(DRHM) for LOCA Licensing Analysis[D].Shanghai:Shanghai Jiaotong University,2011.