XAPR中破口失水事故下堆芯自然循环冷却能力的功能可靠性研究

王宝生，唐秀欢，朱 磊，包利红

(西北核技术研究所，陕西 西安 710024)

西安脉冲堆(XAPR)作为池式研究堆，具有停堆后的非能动安全性，即当反应堆发生事故时，依靠堆池水的自然循环冷却将堆芯余热非能动导出。然而这样的B型非能动系统[1]在缓解严重事故时易受到不确定性的影响，从而引起系统的功能失效[2]。因此，有必要对XAPR堆芯自然循环冷却能力的功能可靠性进行研究，发展一种评价可靠性的方法，为XAPR PSA提供重要的信息。

近年来，Marques等[3]和Jafari等[4]分别提出了功能可靠性评估方法RMPS和REPAS，并用于两相流自然循环系统。谢国锋等[5-6]分别采用响应面(RS)法和蒙特卡罗抽样(MCS)法计算了高温气冷堆余热排出系统功能失效概率。夏少雄等[7]和潘晓磊等[8]分别采用神经网络法和响应面拟合法计算了中国铅基研究堆非能动系统的可靠性。研究表明，功能可靠性评估可采用MCS、减小方差的MCS及RS法。然而上述方法对于小功能失效概率需进行大量的抽样，计算效率不高[9]。

本文结合中破口失水事故(MLOCA)对XAPR堆芯自然循环冷却能力进行可靠性评价。通过改进响应面(IRS)法获得系统输出与输入变量的近似解析表达式，利用重要抽样子集模拟(ISSS)法获得系统输出响应的累积分布函数和置信区间，进而对系统进行功能可靠性评估与分析。

1 功能可靠性评估数学模型

假设影响系统功能的随机输入向量为X=(x1,x2,…,xn)，系统的功能失效概率Pf表示为：

(1)

式中：n为不确定参数个数；g(X)为功能函数，g(X)=A-y(X)，A为失效准则，y(X)表示热工水力学程序，g(X)<0表示系统处于安全状态，g(X)>0表示系统处于失效状态，g(X)=0表示系统处于临界状态；fX(X)为变量X的联合概率密度函数。

为便于计算，式(1)可改写为如下形式：

(2)

式中：I[g(X)]为指示性函数，若g(X)<0，则I[g(X)]=1，若g(X)≥0，则I[g(X)]=0；Ω为整个积分区域。

对于复杂的自然循环系统，难以用快速积分的方法求解式(2)。因为g(X)是一个复杂的热工水力学程序，其功能函数都是模拟程序，需用数值方法求解。由于功能失效概率较小，采用MCS计算需进行多次抽样，且每次抽样都需运行热工水力学程序，计算效率较低。为此，本文提出提高计算效率的两种方法：一是建立响应面替代回归模型，以加快运算速度；二是使用高效蒙特卡罗方法即ISSS法以较少的运算次数获得更加稳健的估计结果。

2 功能可靠性评估分析方法

2.1 IRS法

XAPR堆芯自然循环过程由热工水力学程序模拟，其功能函数均为隐式，可以构建失效响应面的方法来近似表达隐式功能函数[5]。若用响应面法直接求解，需运行热工水力学程序k(2n+1)次，当迭代次数k较大时，计算耗时[8]。为减少计算次数，本文提出IRS法构造失效响应面功能函数，并用该函数代替热工水力学程序。IRS法的流程如图1所示，具体步骤如下。

1) 获得样本点

本文样本点分为两部分：一是通过直接MCS在不确定性输入参数区间内随机抽取M组输入向量，代入热工水力学程序计算相应输出向量；二是在失效域附近进行重要密度函数抽取N组条件样本点，并由热工水力学程序计算相应输出量。

2) 估计响应面功能函数系数

采用bootstrap方法从步骤1获得的样本点中产生t组bootstrap自举样本作为输入样本，并由热工水力学程序计算相应输出量。利用最小二乘法计算出每组响应面功能函数系数，并计算均方根误差，以均方根误差最小为目标函数，确定响应面功能函数系数。

图1 IRS法流程Fig.1 Flowchart of IRS method

3) 功能失效概率计算

从不确定性输入参数的分布中进行ISSS抽样，代入步骤2的响应面功能函数中，获得系统输出响应的累积分布函数和置信区间，进行统计计算得到系统的功能失效概率。

2.2 ISSS法

通过引入合理的中间失效事件，将小失效概率表达为一系列较大的条件概率的乘积形式[10-11]，而较大条件失效概率通过重要抽样来高效估计。给定功能函数g(X)的失效域F，引入一系列中间失效事件临界值0

=…

(3)

式(3)中P(F1)可由MCS计算，条件失效概率P(Fi|Fi-1)可由重要抽样密度函数通过自动分层计算。ISSS法抽样流程如图2所示，具体过程如下。

a——直接蒙特卡罗模拟；b——A1的选取；c——重要抽样模拟；d——A2的选取图2 ISSS法抽样流程Fig.2 Illustration of importance sampling subset simulation procedure

i=1,2,…,m-1

(4)

5) 重复步骤3和4，直到满足条件Am≥A，此时自动分层结束。

6) 计算功能失效概率估算值为：

(5)

3 XAPR堆芯自然循环功能可靠性分析

3.1 系统安全功能

XAPR采用池式反应堆结构，堆芯采用自然循环冷却。当反应堆停堆后，堆芯余热可依靠池水的自然循环冷却，池水与环境间的热量交换即可将堆芯余热非能动导出。当发生失水事故时，事故进程可划分为3个阶段：1) 水冷自然循环冷却阶段；2) 堆芯半裸露阶段；3) 空气冷却自然循环阶段。这样即可在失水事故时建立堆芯长期自然循环冷却，将堆芯余热导出，保证燃料元件的完整性。

3.2 系统建模与验证

采用RELAP5模拟失水事故堆芯自然循环冷却过程。XAPR热工水力建模计算节点示于图3。图3中控制体101～107为堆池，120～124为堆池混凝土壁外空气边界，150～163为一回路管道，170～172为简化二回路边界，109为堆池顶部空气边界。为验证模型的准确性，采用稳态运行工况作为边界条件进行计算，得到堆芯自然循环流量为12.21 kg/s，该流量与实际自然循环流量12.13 kg/s基本一致。

3.3 系统不确定性参数

XAPR自然循环运行时涉及到2种不确定性：第1种是“偶然性”所导致的不确定性，它与模型的几何性质有关；第2种是“认知局限”所导致的不确定性，它与运行和实验数据的缺乏而导致的对有关现象认知的局限性有关[10]。本文主要研究因“认知局限”所导致的不确定性。

图3 RELAP5计算节点图Fig.3 Nodalization scheme of RELAP5 calculation

通过层次分析法选出14个不确定性输入参数。在目前处理中，当可得到的数据有限时，各输入参数的概率分布及取值区间主要以工程设计标准为基础，并结合专家评价的主观方法得出。表1列出输入参数的概率分布及特征参数。这些参数的取值及分布考虑了所有可能的事故工况，是一种较为保守的方法。

3.4 系统失效准则

定义功能准则为：当包壳温度低于或等于500 ℃时，燃料芯体最高温度应低于1 150 ℃；包壳温度高于500 ℃时，燃料芯体最高温度应低于970 ℃。从概率安全角度出发，假设燃料芯体最高温度超过安全限值970 ℃则认为功能失效。设向量X为系统不确定性参数(如表1所列)，To,max(X)为燃料芯体最高温度，则功能函数g(X)可表示为g(X)=970-To,max(X)，当g(X)<0时，即认为失效。

4 功能可靠性评价与结果分析

4.1 燃料芯体最高温度置信水平

用于拟合响应面功能函数的样本点包括：第1部分采用MCS抽取M=120组随机样本点；第2部分采用重要密度函数抽取N=24组样本点；第3部分在接近失效域附件采用bootstrap方法抽取12组样本点。共156组样本点，其RELAP5运行计算耗时为159.8 h，其中tM=120.1 h、tN=27.3 h及ttb=12.4 h。利用156组样本点求解响应面功能函数的系数，结果如下：

1.72Tg+0.18p-31.19ξ1+1.34×102ξ2-

6.3ξ3+6.21ξ4+0.09dx+13.48kin+

12.98kout-4.79td+8.54tm+12.89tp-

(6)

图4示出燃料芯体最高温度对比与残差散点图。由图4a可看出，响应面函数的拟合值与REALP5的输出值吻合很好。由图4b可看出，各工况的残差均在±6 ℃以内。因此可将响应面函数作为RELAP5替代模型进行功能失效概率计算。

为研究IRS在抽样效率上的提升，计算了156组抽样在RS和IRS抽样中的比例。在利用156组抽样得到的输出结果中，RS大约有2组数据导致To,max大于970 ℃，所占比例约为1.28%；IRS得到15组数据导致To,max大于970 ℃，所占比例约为9.61%，失效事件出现的概率提高了将近8倍，抽样效率明显提升。图5示出RS和IRS对20组样本点的To,max的预测结果。

表1 输入参数的概率分布及特征参数Table 1 Probability distribution and characteristic parameter of input parameter

图4 燃料芯体最高温度对比(a)与残差散点图(b)Fig.4 Comparison of fuel peak temperature (a) and residual scatter plot (b)

图5 20组样本点的To,max预测结果Fig.5 Prediction result of To,max of 20 samples

图6示出α=0.9、0.95、0.99及0.999等4种不同置信水平下燃料芯体最高温度的累积分布函数(CDF)。由图6可看出，不同置信水平下，ISSS均有相当数量的样本点落在计算区域内。但MCS抽取400个样本点时，在高置信水平下，很难得到精确值。这是由于当抽样数为400时，只有平均400×(1-0.99)=4个样本点落在0.99～1区间内，如图6c所示；而只有400×(1-0.999)=0.4个点落入0.999～1区域，如图6d所示。然而，对于ISSS能保证有100、10个样本点落入相关区域，这就给可靠性分析提供了足够的信息。注意到MCS抽样10 000个样本点，分别也有100、10个样本点落入目标区域，如图6c、d所示，但抽样次数较多，计算效率较低。

4种不同置信水平下，燃料芯体最高温度均值与标准差列于表2。由表2可看出，ISSS的结果与MCS抽样10 000次的结果基本一致，即便是在高置信水平下，两者的标准差也相差不大，但ISSS的抽样数明显低于MCS的抽样数，仅约为后者的1/25。

α：a——0.9；b——0.95；c——0.99；d——0.999图6 不同置信水平下燃料芯体最高温度的累积分布函数Fig.6 CDF for fuel peak temperature at different confidence levels

表2 不同置信水平下燃料芯体最高温度均值与标准差Table 2 Mean value and standard deviation of fuel peak temperature at different confidence levels

4.2 功能失效概率计算及比对

不同层次下燃料芯体最高温度的自动分层结果如图7所示。由图7可看出，各中间失效事件F1、F2、F3、F4的临界值A1=843.52 ℃、A2=910.25 ℃、A3=948.79 ℃及A4=999.08 ℃，由于A4=999.08 ℃>970 ℃(图7d)，自动分层结束。根据式(5)计算自然循环冷却能力的功能失效概率为Pf=3.796×10-3。为便于各种计算方法的比较，还同时采用了MCS、IRS、重要抽样(IS)法进行计算，计算结果列于表3。

由表3可知，ISSS+IRS抽样400次，即可达到MCS+RELAP5抽样104时的计算精度。MCS+IRS抽样400次时，其计算结果为1.875×10-3，导致这样的结果是因为当失效概率为10-3数量级时，只有约400×10-3=0.4个样本点落入失效域内。IRS实质是曲线拟合近似方法，其计算精度依赖于功能函数的线性程度，XAPR自然循环的功能函数非线性程度较高，所以结果欠佳。MCS+RELAP5抽样效率低、计算费时，在抽样次数为104的情况下计算耗时约为ISSS+IRS的6.25倍，因而不适用。在相同抽样次数的情况下，IS+IRS计算耗时与ISSS+IRS接近，但其计算精确度略差。

4.3 参数灵敏度分析

灵敏度分析可反映输入参数的不确定性对系统功能可靠性的影响，从而为提高系统功能可靠性提供指导性建议，有效地减小不确定性。由于非能动物理过程为隐式非线性非单调关系，传统表征参数灵敏度的相关系数(CCS)、标准回归系数(SRCS)、秩相关系数(RCCS)及标准秩回归系数(SRRCS)等具有局限性[12]。本文将基本变量xi(i=1,2,…,n)均值μxi的变化引起的功能失效概率变化的比率表征为灵敏度，在数学上表达为由失效概率Pf对基本变量均值μxi的偏导数，即∂Pf/∂μxi。为消除量纲影响，引入了无量纲灵敏度系数Sμxi=(∂Pf/∂μxi)×(σxi/Pf)[13]，其中σxi为输入变量标准差。灵敏度分析结果如图8所示。

a——层次1；b——层次2；c——层次3；d——层次4图7 不同层次下的燃料芯体最高温度Fig.7 Fuel peak temperature at different conditional levels

表3 不同计算方法计算的功能失效概率Table 3 Functional failure probability for different calculation methods

由图8可见，空气入口温度Tg、衰变热功率相对偏差Q、流道进口阻力系数kin、堆池水入口温度Tw及排水泵开启延迟时间tp对系统输出有较大影响，其中Tg、Q及kin最为敏感。因为这些输入参数直接影响事故发生后堆芯自然循环的输热能力和建立时间，对堆芯余热排出功能可靠性影响较大，因此，减小上述输入参数的不确定性可有效降低堆芯自然循环的功能失效概率，提高系统功能的可靠性。

图8 输入参数的无量纲灵敏度系数Fig.8 Normalized sensitivity coefficient of input parameter

5 结论

1) 本文方法将IRS和ISSS相结合，基于失效响应面函数将功能失效概率表述为一系列较大条件失效概率的乘积，通过子集模拟使得失效域内样本计算效率大为提高。相比于其他方法，该方法具有较高的抽样效率，同时又能保证较好的计算精度，且对于非线性程度较高的自然循环系统有很强的适应性。

2) XAPR中破口失水事故时，由“认知局限”不确定性导致的堆芯自然循环失效总有非零的发生概率，其失效概率为3.796×10-3。为此，在涉及非能动自然循环的可靠性评估中，功能失效是重要的，其失效概率应基于热工水力计算和可靠性方法进行全面评价。

3) 减小Tg、Q、kin、Tw及tp等5个关键参数的不确定性，可更有效地降低XAPR堆芯自然循环的失效概率，提高XAPR系统的功能可靠性。