曹志伟,林支康,王 婷,梁 任,鲍 杰,卢向晖
(1.中广核研究院有限公司,广东 深圳 518031;2.生态环境部 核与辐射安全中心,北京 100082)
大破口失水事故(LB LOCA, large break loss of coolant accident)是核电厂最重要的设计基准事故之一,是验证应急堆芯冷却系统(ECCS)能力并直接限制核电厂发电能力的极限事故,LB LOCA裕量也是最重要的核电厂安全裕量指标之一。由于LB LOCA现象复杂,后果严重,在核电厂安全分析中处于非常重要的地位。LB LOCA分析程序、方法及分析均是核电厂建造和装料许可证等申请的关键内容。目前拥有的LB LOCA分析方法是20世纪90年代开发的较为保守的分析方法,保守的分析方法掩盖了核电厂的真实安全裕量,如18个月换料模式下的CPR1000核电厂,最低LOCA裕量只有1%左右。保守分析方法无法真实反映核电厂设计的安全指标和先进性,也极大地限制了核电厂经济性的提升。因此,为提升核电厂安全裕量,本文基于现有的最佳估算热工水力分析程序CATHARE GB及相关保守方法,开发基于统计理论的先进LOCA分析方法,并初步应用于CPR1000核电厂。
20世纪90年代,法国电力公司和法玛通公司基于其多年的核电厂设计和运营经验,将LB LOCA分析方法逐步从保守的、满足美国联邦法规10CFR50附录K[1]的方法论转向了基于最佳估算程序的分析方法——确定论现实方法(DRM, deterministic realistic method)[2]。程序计算结果的不确定性来源于两方面:初始边界条件和软件物理模型的不确定性。分析所用的计算程序是现实的,即程序必须对所有瞬态物理现象进行最佳估算。DRM开发基于CATHARE 2 V1.3L程序,针对LB LOCA现象提供了专用的模型。这些模型均与现有实验进行对比,得到程序计算模型的不确定性。基于现实的电厂模型,通过对程序计算模型的不确定性、电厂技术参数及初始状态的不确定性进行抽样计算,得到每个目标参数(如燃料包壳峰值温度,PCT)在95%置信水平下、95%概率下不会被目标参数值超过的保守值(即双95%值)。随后定义出程序的保守模型(DRM惩罚模型),该模型计算结果可包络目标参数的双95%值,因此通过DRM得到的目标参数结果高于统计方法得到的95%置信水平的结果。所构造的惩罚模型不仅要求保守性最小,还应保留程序对真实瞬态工况的模拟能力。DRM的开发主要包括4个方面:1) 评估模型的现实性;2) 量化不确定性;3) 引入惩罚模型;4) 保守性评价。
在DRM中为了满足计算结果的保守性,引入了程序惩罚模型。为了尽可能地既不影响瞬态过程,又能同时包络PCT95/95值,最终选择对几个参数进行惩罚,称为DRM惩罚模型[1]。
在LOCA保守评价模型中,影响LOCA评估的相关模型和电厂状态依据美国联邦法规10CFR50附录K[1]的要求均须采用保守方式处理。在LOCA最佳估算评价模型[3-4]中,程序物理模型和电厂状态参数的不确定性采用统计方法处理,以综合呈现出主要评估参数的不确定性范围。影响LOCA最佳估算不确定性的因素分为两大类:一类是程序模型的不确定性,另一类是电厂状态参数的不确定性。
中广核确定论统计方法(GSM, CGN deterministic statistic methodology)为一种介于保守评价模型和最佳估算评价模型两类方法之间的折衷式LOCA分析方法。在该方法中,重要程序模型部分以经实验验证的保守方法处理;而电厂重要状态参数的处理,则遵循CSAU的执行程序,采用统计方法量化电厂重要状态参数不确定性的整体效应。对于无法实行不确定性分析的电厂重要参数则采用包络值,以保守处理该参数的不确定性影响。
GSM是基于现有的最佳估算热工水力分析程序CATHARE GB以及DRM分析相关惩罚模型假设而开发的LOCA分析方法,在GSM中,考虑的不确定性主要有3类:1) CATHARE GB程序模型相关不确定性;2) 电厂模型相关假设;3) 电厂状态相关参数的不确定性。
对于CATHARE GB程序模型相关不确定性,采用了DRM惩罚模型,相关参数和取值列于表1。表1中:HTCBerenson为喷放阶段堆芯传热系数;HTCReflooding为再淹没阶段堆芯传热系数;Tmsfb为最小膜态沸腾转变温度;Qle为气液两相间的冷凝系数。
表1 惩罚模型相关参数和取值Table 1 Relevant parameter and value of penalization model
电厂模型相关的不确定性包括:考虑单一故障准则;破口发生时刻考虑丧失厂外电源,主泵停运;事故发生时的燃料循环,采用保守的燃耗数据,通过敏感性分析得到;一回路每个环路的流量为热工水力设计流量;考虑最大的堆芯旁流量,使通过堆芯的冷却剂量最少;一回路质量计算考虑蒸汽发生器堵管率;低安全壳压力使得再淹没进程减缓并降低再淹没期间堆芯内热传递系数,对PCT更加不利,因此在计算中采用正常运行工况下最小的初始安全壳温度和压力,并选择安全壳最大湿度。破口尺寸和类型采用统计抽样处理,其中,破口类型为等概率分布,破口面积的抽样范围为[0.241 6A,2.0A],A为冷管段的横截面积。
电厂状态参数的不确定性处理重点是要筛选出重要电厂状态参数并对其进行不确定性量化处理。
GSM对于电厂状态相关参数采用了统计抽样处理,针对CPR1000核电厂,采用抽样统计处理的电厂状态相关参数列于表2。
表2 重要电厂参数及其不确定性范围和分布形式Table 2 Plant status parameter and it’s uncertainty range and distribution
用统计方法量化计算误差时,必须针对相关不确定性因素进行随机取样,产生足够数目的随机样本。随机样本数目由所采用的统计方法决定。一般而言,统计可分为非参数方法和参数方法两大类。非参数方法的统计分析不需要知道分析结果误差概率密度分布类型,而参数方法的统计分析必须假设分析结果误差概率密度分布特性,如函数种类、均值及标准方差等[5]。GSM采用简单随机抽样,对电厂参数进行统一抽样,得到59组统计抽样组合工况并计算得到59个PCT值,供上述两种统计方法分析使用,计算PCT95/95统计上限值。
采用考虑DRM保守惩罚的CATHARE GB程序模型对59组工况计算得到59组PCT,并对此进行统计处理以量化在考虑电厂状态参数不确定性后的整体计算结果的误差范围,并进一步计算或估算出PCT95/95统计上限值。采用的统计处理方法有两种,即正态分布单侧置信上限方法(参数统计法)和顺序统计方法(非参数统计法)。
1) 正态分布单侧置信上限方法[6]
在此统计分析方法中先假设59组工况的分析结果满足正态分布,并经适当方法检验确认,如χ2拟合检验法。如果满足特定函数分布,则可由样本均值μs及样本标准方差σs在一定的置信度95%下估算总体均值μp及标准方差σp,并进一步估算上限值,如PCT95/95之值。样本均值和样本标准方差分别为:
其中:xi为样本变量参数值;n为样本数量。
参数统计分析可分为3个步骤,即正态分布检验、由样本均值与标准方差估算总体均值与标准方差及由分析结果计算95/95统计上限值。
(1) 正态分布检验
所采用的正态分布检验方法是χ2拟合检验法,该方法有两项适用要求:只适合大样本情形,一般要求样本容量N≥50;每个区间的样本数(NPi)不应太小,其中Pi为子区间i对应的概率,一般要求不小于5,否则应适当合并邻近区间,使得NPi≥5,同时区间数k一般应在4~20之间。当假设分布为正态分布N(μ,σ)时,真实的总体μ、σ仍是未知。根据极大似然估计法可由样本μs和σs合理估算。
将样本数NPi依大小分为数个区间,定义正态分布的分歧度υ,量化各区间预期样本数和实际样本数的差异:
(2) 由样本均值与样本标准方差估算总体均值与标准方差
以PCT上限包络值的计算为例,总体均值μp在1-α置信度下单侧置信上限为:
Pμp≤μs+tα(n-1)σs/n{}=1-α
其中,tα为学生氏分布统计参数。
同样地,总体标准方差σp在1-α置信度下单侧置信上限为:
(3) 计算95/95统计上限值
上述相关计算式中在一定置信水平或置信度1-α下,μp和σp上限包络值可由样本均值和标准方差直接计算。当μp和σp上限包络值在95%置信度下计算出后,目标参数的95/95统计上限值Y95/95可直接由下式计算:
Y95/95=μp,95%+1.645σp,95%
其中,μp,95%和σp,95%分别为在95%置信度下μp和σp的上限包络值。
2) 非参数统计法
采用Wilk’s非参数统计法[7]不需先行确定样本所服从的概率分布,而是在有限的随机抽样样本中,如n=59、93或124组,保守估算可能的95/95统计上限值:Y95/95≈Y1st(59)、Y95/95≈Y2nd(93)或Y95/95≈Y3rd(124)。
3) 最终目标参数值的确定
在确定最终PCT值时,如果正态分布单侧置信上限方法和非参数统计法均能计算或估算出PCT95/95值,用较保守者作为最终PCT值,即:
PCT=max(PCT95/95,PCT非参数统计)
将GSM应用于CPR1000核电厂某燃料管理方案研究项目中,其分析结果与传统保守的DRM进行比较,以合理挖掘LB LOCA的裕量。不确定性分析的重要电厂参数范围和分布列于表2。经统一随机抽样59组样本组合工况,并使用CATHARE GB程序模型计算得到59组样本对应的PCT值。GSM应用的流程如图1所示。
根据Wilk’s公式γ=1-βn,当概率水平(γ)、置信水平(β)均设定为95%时,n即为59。此时59组抽样结果的最大值即为95/95统计上限值的保守估算结果。采用Wilk’s非参数统计方法计算得到PCT95/95统计上限值为1 136.63 ℃。
参数统计分析分为3个步骤:1) 正态分布检验;2) 由样本均值与标准方差估算总体均值与标准方差;3) 由分析结果计算95/95统计上限值。
1) 正态分布检验
图1 GSM流程Fig.1 Procedure of GSM
首先样本均值和标准方差的计算如下:
2) 由μs和σs估算总体均值与标准方差
μp的95%单侧置信区间为(-∞,899.244],σp的95%单侧置信区间为(-∞,94.514]。
由于满足正态分布,所以95/95的PCT值为:
PCT95/95=μp,95%+1.645σp,95%=
899.244+1.645×94.514=1 054.72 ℃
表3 拟合优度检验分歧度表Table 3 Data for goodness-of-fit test
3) 计算95/95统计上限值
由Wilk’s非参数统计法得到PCT1st=1 136.63 ℃,最终计算得到PCT95/95的值为:
PCTGSM=max(PCT1st,PCT95/95)=1 136.63 ℃
GSM应用于CPR1000核电厂LB LOCA分析结果与传统DRM分析结果的比较如图2所示。由图2可见,采用GSM计算得到的LB LOCA分析结果比传统DRM分析结果低117 ℃,可提高约9%的裕量。
图2 GSM与DRM分析结果对比Fig.2 Comparison of analysis result of GSM and DRM
GSM是介于保守评价模型和最佳估算评价模型两类方法之间的折衷式LOCA分析方法。在该方法中,程序物理模型仍以DRM惩罚模型进行保守方法处理,对核电厂重要状态参数采用统计方法量化各重要电厂状态参数不确定性范围和分布,并对统计抽样计算得到的目标参数(PCT)进行统计处理以得到PCT95/95。GSM应用于CPR1000核电厂LB LOCA分析可挖掘约9%的LOCA裕量。