杜 芸,李焕鑫,梁国兴
(上海交通大学 核能科学与工程学院,上海 200240)
风险指引安全裕度特性分析(Risk Informed Safety Margin Characterization,RISMC)方法论旨在开发一套先进的风险评估方法,耦合风险指引型的安全裕度的量化,最终能够用于核电厂的决策制定[1-3]。为获得风险指引安全裕度支持决策者制定决策计划,核电厂安全分析方法不能再局限于分析设计基准事故序列,应该对所有的概率显著的序列进行考虑和分析。下一代安全分析技术需要具备的核心技术和能力主要包括[4]:(1)事故序列的生成、事故序列的特性描述和量化;(2)概率显著序列(probabilistically significant sequences,PSSs)的现实性现象分析。这提示我们先进的安全分析方法需要对原有的两种分析方法——确定论分析方法(DSA)和概率论分析方法(PSA)[5]进行整合和改进(见图1)。
图1 以RISMC方法论为指导的新一代分析技术Fig.1 Next generation analysis capability based on RISMC
风险指引的安全分析方法论可提供三类主要的应用,其分別为:(1)量化风险指引的安全裕度,(2)从风险的角度优化核能电厂各重要系统的设计,以及(3)量化系统老化或设计变更对系统所造成的风险响应。本文将以RISMC方法论为指导开发一套切实可行的先进安全分析方法,主要关注的是第三类应用中的小幅功率提升对电厂全厂断电事故风险的影响。对于小幅度的功率提升,传统的PSA方法由于在成功准则的设定上过于保守并且没有综合考量认知性参数和随机性参数的不确定性等缺陷,所以无法量化其风险响应的变化。但是任何的设计变更对电厂造成的风险冲击都必须被量化,所以该方法的开发十分的有必要。
本文首先基于RISMC方法论开发一套全新的安全分析方法流程,然后利用该方法对真实电厂的全厂断电事故的重要序列进行堆芯熔毁概率(CDF)的重新计算和评估,并和传统PSA的分析结果进行比较,然后进一步的利用该方法量化功率小幅提升对电厂风险的影响。
将传统的PSA方法和最佳估算加不确定性方法(BEPU)[6]整合起来,即保留DSA中对事故现象的描述功能和对相对全面的不确定性进行量化的功能,同时也发扬PSA中对事故多样化序列与事故场景进行风险量化的功能。通过传统PSA方法的分析,所有对目标始发事件中的概率显著序列即风险最显著的序列将会被挑选出来,并且建立改进的PSA模型;利用BEPU技术,随机性的不确定性参数以及认知不足性的不确定性参数将会通过对重要现象的分析研究而被挑选出来,并且被量化。随机性的参数主要包括系统或元件的随机失效组合,外电恢复的时间,系统或元件失效的时间或者操作员决策时间等;认知不足性的参数主要包括电厂的运行参数,模型的不确定性参数等。
当使用了对两类参数进行同时随机抽样的技术同时松开对成功准则的设置之后,堆芯熔毁概率将会是三个部分的乘积,分别是始发事件的发生概率,序列的发生频率以及条件失效概率,如公式(1)所示。结合了确定论不确定性分析思想与技术的PSA模型将会更加具有现实性。条件失效概率将会是由概率分布得到的一个比值,如此计算得到的安全裕度将带有概率的属性,具有风险指引的意义,对电厂决策的制定更加有利。
CDF=Fie·Pseq·Pce
(1)
由于各类参数不确定性抽样之后,需要在风险评估过程中执行大量的系统程序的模拟计算,以量化参数不确定性的统计结果,所以将这种基于风险指引的安全裕度特性分析方法论的先进安全分析方法取名为计算风险评估(Computational Risk Assessment,CRA)方法。
由于下一代的安全分析方法要求能够高效地完成风险量化。而对分析过程中的抽样范围进行合理的缩小,不仅能够提高分析的敏感性,更能提高分析的效率。尤其对于随机性的不确定性参数,比如外电恢复的时间或者操作员做决策的时间等。比如,在SBO传统PSA模型中的关键前沿事件“外电是否及时恢复”就可以被划分成三个部分,如图2所示,即当外电在tmin之前恢复的情况(确定成功),外电在tmax之后回来的情况(确定失败)以及外电在tmin和tmax之间恢复(条件失败)三种情况进行讨论。其中,tmax定义为当所有电厂状态参数都在最良好的状态时 (best bounding case),满足堆芯不熔毁条件的外电恢复的时间点;tmin定义为当所有电厂状态参数都处在最不理想的情况时(worst bounding case),满足堆芯不熔毁条件的外电恢复时间点,如图3所示。当外电的恢复时间落在小于tmin和大于tmax这两个时间段时,电厂状态是确定的,没有必要进行外电恢复时间的抽样。而外电恢复时间落在(tmin,tmax)时,需要进行抽样,并且运用系统程序计算得到真实的电厂状态。这样可以有效缩窄抽样的区间以及强化抽样个案的有效响应,亦可节省抽样的数量以及程序计算的时间,将抽样样本的数量能够高效地运用于电厂运行结果不确定的tmin和tmax之间的时间段上。
图2 电力恢复时间的三个时间区间Fig.2 Time intervals of three subheadings in power recovery time distribution
图3 丧失热阱为主因序列的tmin和tmaxFig.3 The tmin and tmax in loss of heat sink dominant group
综合以上分析,开发以RISMC方法论为指导的计算风险评估(CRA)方法,具体的实施步骤如下所示。图4对CRA方法的流程进行了描述。
图4 计算风险评估方法的流程图Fig.4 Process of the computational risk assessment
(1)确认目标始发事件,量化其发生概率
对于特定的设计变更,只有特定的始发事件会受到其影响。
(2)确认始发事件下的概率显著序列
对特定始发事件进行PSA模型的计算分析,得到概率风险显著的序列(PSSs),挑选出需要进行细致分析的概率显著序列。
(3)建立改进的PSA模型并且量化序列的发生概率
为了对选定的概率显著序列建立一个改进的PSA模型,使得该改进的PSA模型可以更高效敏感的处理随机参数的不确定性,需要对涉及随机参数的前沿事件进行进一步的细分。其中需要考虑实际情况以及运用预测的极限面[7](limit surface)。这种前沿事件的细分对于随机参数来说,也是一种数值上的划分,能够缩小抽样范围。
(4)针对概率显著序列,确定主要的随机性参数和机械性参数
针对改进的PSA模型中的概率显著序列进行事故序列的过程分析,从中得到对结果有明显影响的随机性不确定性参数以及机械性不确定性参数。同时探究这些参数不确定性的来源,参数的分布特点以及分布范围。最终确定不确定性参数的抽样范围以及抽样方法。
(5)对随机性参数以及认知不足性参数进行随机抽样,组成N个抽样个案
随机性参数以及认知不足性参数都需要使用蒙特卡罗方法[8]或者拉丁超立方方法[9]进行随机抽样,每个参数之间都是相互独立的。每抽取一组参数就形成一个试算案例,直至抽取到满足抽样要求的指定采样个数N,即组成N组试算案例,用来模拟概率显著序列的事故发展。
(6)运用合适的系统程序对每一组抽样形成的案例进行模拟计算
为了得到条件失效概率,每一组试算案例的PCT值都可以通过合适的系统程序计算得到,例如RELAP5/MOD3。经过N次的程序计算,最终所有N个PCT的值将会被计算得到。
(7)得到N个PCT所对应的概率分布,从而得到PSS对应的条件失效概率
由于用代数方法求得条件失效概率的值往往被抽样样本数所左右,其收敛性还需要进一步探究和证实。在考虑效率以及抽样个数不宜过大的同时,决定采用统计学的方法检验或者拟合得到N个PCT值所服从的概率分布函数,而后得到指定PSS的条件失效概率。
(8)对每一个PSS进行堆芯熔毁概率的计算
利用步骤1中获得的始发事件发生概率,步骤3中获得的事故序列发生的概率以及步骤7中获得的失效概率相乘得到目标概率显著序列的堆芯熔毁概率。
针对小幅功率提升这一特定设计变更,其对SBO事故会造成风险的冲击,所以对SBO进行风险量化分析。由传统的PSA模型[10]给出的计算结果,某传统三环路压水堆的SBO事故中的概率显著序列是序列A和序列B(见图5)。其中序列A为电厂发生SBO事故后,主泵轴封没有发生早期失效情况下电力没有及时恢复导致的堆芯熔毁序列;而序列B为轴封发生早期失效叠加电力没有及时恢复的序列。本文运用新的安全分析方法评估序列A的堆芯熔毁概率,今后的工作中会继续对序列B进行重新评估计算,以获得更加完整全面的结果。根据CRA方法的指引,对传统的PSA模型进行改进。
图5 含有计算结果的全厂断电事件树Fig.5 Event tree for station blackout with probability
针对传统SBO模型中的序列A和序列B,改进主要体现在三个方面:第一个是序列A和序列B分界点的重新界定,即判断轴封早期失效的时间点;第二个是关键前沿事件的细分;第三个是引入抽样技术,计算条件失效概率。
对于一个典型压水堆的全厂断电事故序列A,反应堆立即成功停堆。初期,反应堆冷却系统的边界是完整。汽动辅助给水系统能够成功的运行8个小时。由于缺少了所有的交流电源,冷却泵将会失去轴封冷却水的注入,会在全厂断电发生后立即发生轴封漏水。每个冷却泵的轴封漏水流量为21 GPM。进一步的,由于冷却泵长期缺少轴封冷却水的冷却,其会发生轴封失效,从而每个泵会产生流量高达450 GPM的冷却剂流失[11]。三个冷却泵轴封失效所产生的冷却剂流失相当于冷管段直径为一英尺的破口所产生的流量,这相当于冷管段的小破口流量。在传统PSA模型中,轴封早期失效的判断时间点是由非常简化和守的计算得到,在改进的PSA模型中,采用系统程序随电厂建模(见图6)分析得到。额定工况下,全厂断电且轴封发生泄漏但是不发生失效的情况下,为避免堆芯熔毁安注系统最晚的投入时间为14.23 h。该时间将作为轴封是否早期失效的判断时间。所以在新模型中,失去热阱为主因的试算案例(轴封失效时间大于14.23 h)得到的堆芯熔毁概率将与传统PSA中的序列A的堆芯熔毁概率(2.96×10-7)进行对比。
对传统的PSA模型(见图6)的关键前沿事件进行细分,这样的细分会增加整个方法的高效性,将抽样和计算集中到最有效的范围内。由改进PSA模型图(见图7)中可以看到,序列A堆芯熔毁的概率对应新模型中的序列3和序列5的堆芯熔毁概率之和。
而电源恢复的时间确定之后,由于电厂状态存在不确定性,电厂安全到底是何种走向必须由系统程序的模拟计算结果决定。实际情况下,必须考虑电厂状态的不确定性以及随机性参数的不确定性,所以最晚的外电恢复时间无法用一个单一的准则来确定。在改进的PSA模型中,对认知性参数和随机性参数抽样将会决定多组试算个案,每一个个案代表着一种电厂有可能的工况,运用系统程序对每一种工况进行单独的模拟计算,得到每个工况下的电厂终态。最终,针对N组算例的结果,序列5的条件失效概率需要采用统计的手段得到,不再是0或者1。
图6 典型三环路压水堆建模节点图[12]Fig.6 Nodalization of a typical three-loop nuclear power plant
图7 基于计算风险评估方法的改进的PSA模型事件树Fig.7 Revised event tree based on computational risk assessment
经过名义工况下对全厂断电事故的模拟分析可知,对于最终包壳峰值温度来说,首先重要的是对象的初始状态,即电厂进入SBO时的状态。其中包括内热源的功率,堆芯初始的温度,汽动辅助给水的流量温度,蒸汽发生器的水位以及电厂内的压力等。其次重要的是在电力恢复之后,电动辅助给水系统和安注系统的启动是否能将堆芯内的温度迅速降下来,以防堆芯发生熔毁,此时必须关注给水以及安注水的流量和温度等。除此之外,在全厂断电事故中有两个最为重要的时间点。一个是轴封失效的时间,因为轴封的失效对于反应堆来说相当于经历一个小破口事故,这对于事故进程来说是十分不利的。轴封失效的时间根据传统PSA中的数据,满足威布尔分布。另一个重要的时间点是外电(包括柴油发电机)恢复的时间,电力恢复的时间通常满足对数正态分布。本文研究对象电厂的电力恢复时间通过贝叶斯方法处理满足如下对数正态分布。
综合对事故过程的分析以及专家的意见,最终确定对8个电厂状态参数,以及3个随机参数,具体参数如表1,表2所示。
表1 电厂状态参数的抽样范围和分布
表2 随机参数抽样范围和分布函数
为了让改进的PSA模型对设计变更更加的敏感,将有效的电力恢复时间的抽样范围需要缩小到tmin和tmax之间,其中tmin和tmax是由最差工况和最佳工况决定的。对于失去热阱为主因的事故序列A,最佳工况和最差工况的取值如表3所示。
表3 丧失热阱为主因序列的最佳工况和最差工况取值表
由系统程序对两种工况下全厂断电并且电力及时恢复的过程进行模拟和分析(见图3)。对应失去热阱为主因的序列A,tmin和tmax分别为13.071 h和13.447 h,因此(13.071~13.447)即为外电恢复时间抽样范围;两种电源的恢复时间需要同时抽样,任何一种电源的恢复均可支撑电厂的救援系统启动和运行。
根据对传统PSA模型的改进,对相关随机性参数和机械参数的选择以及其分布和抽样范围的确定,一个基于CRA方法的新的PSA模型被建立(见图7)。
由公式(1),堆芯熔毁概率的计算需要由三个部分的组成。根据改进的PSA模型以及确定的分支点,序列的发生频率可以通过计算得到。
当轴封失效时间大于14.23 h时,试算案例被划分为失去热阱为主因的事故序列。通过预先确定的临界轴封失效时间和条件成功抽样时间区间(tmin,tmax)以及抽样概率密度分布函数的确定,可以直接计算得到SL、BS、BF和BC分支的概率(见图7)。轴封失效时间大于14.23 h的概率是0.006 9。对于燃气轮机在0时刻启动和成功运行的事件没有发生概率是0.27。
为了与传统PSA模型中的序列A进行对比,在改进PSA中,需要计算序列3、4、5的发生概率。根据图7,除了已经在上文中介绍过的部分以外,还需要知道电力恢复一定能阻止堆芯熔毁事件的互斥事件的发生概率以及在此情况下电力恢复一定不及时事件的概率。对于图7中的“确定成功”,其具体意义为“至少有一种电源的恢复是早于tmin”。该概率通过表2给出的概率分布函数以及公式(2)计算得到为0.984 6。所以互斥事件发生的概率为0.015 4。
P(Bs)=p(tosr
p(tDGr>tmin)
(2)
而电力恢复一定不及时的概率为两种电源的恢复时间均大于tmax,所以Bf事件发生的概率为0.927 8,见公式(3)。由公式(4)得到S3发生的概率为2.662×10-5。
(3)
(4)
由公式(5),计算得到序列4和序列5两个序列发生的概率之和为2.07×10-6。
(5)
改进PSA模型关于S4和S5最后一个分支的概率需要对抽样产生的试算案例的模拟计算结果进行分析评估才能得到。这一点同时也是改进的PSA模型中最为灵活以及符合实际现象的体现。
对10个重要参数进行随机抽样,抽样样本数为50。8个电厂状态参数采用最经典的蒙特卡罗抽样方法。对于3个随机参数,参数的分布函数以及具体情况较复杂,所以采用更适合小样本抽样的拉丁超立方方法进行抽样。将所有抽样得到的参数随机组合,形成50个抽样个案。采用RELAP5/MOD3对50组试算案例进行事故模拟,50条随时间变化的PCT曲线如图8所示,50个PCT值在散点图上表示,如图9所示。
图8 丧失热阱为主因序列算例的50个包壳峰值温度随时间变化曲线图Fig.8 PCT responses of 50 trials of theloss of heat sink dominant group
图9 丧失热阱为主因序列算例的50个包壳峰值温度结果散点图Fig.9 50 sets of PCT results in lossof heat sink dominant case
50个PCT值在数学意义下的平均值和方差可以通过最大似然估计的理论获得:
(6)
(7)
假设H0:该组数据服从正态分布,分布函数为N(1413.18,191.892)。
表4 卡方检验分组和计算结果
其中,pi是基于假设函数落在对应区间的频率;npi为基于假设的区间内的预期个数
=7.815>5.522=χ2
(8)
图10 丧失热阱为主因序列算例的结果直方图Fig.10 50 sets of PCT results inloss of heat sink dominant case
根据正态分布数据的特性[14],在一定的置信度下,总体的平均值和标准差的估计如公式(9)和公式(10)所示。
(9)
(10)
得到当置信度为95%,总体服从正态分布且均值和标准差为(1 458.68,230.6)(见图10)。所以PCT大于1477.6 k的概率pce,5为46.73%。
序列3和序列5的堆芯熔毁概率由公式(11)和公式(12)计算得到为3.162 3×10-8每堆年和1.149×10-9每堆年。
CDF(s3)=Fie·P(s3)·Pce,3
(11)
CDF(s5)=Fie·P(s4&s5)·Pce,5
(12)
CDF(revised SA)=CDF(s3)+CDF(s5)
(13)
综合序列3和序列5的结果,改进下的序列A的堆芯熔毁概率为3.28×10-8每堆年(见公式(13))。对比传统方法的2.96×10-7,CRA方法下的结果缩小约九倍。这是由于传统PSA模型中轴封早期失效的临界时间是按照保守的方式确定的,与本文中利用系统程序计算的方式不一样,导致该分支概率有较大的不同。如果除去该影响因素,CRA方法的计算结果较原来的方法缩小了28% 。当然,在后续的工作中,改进下的序列B的堆芯熔毁概率也将被量化,综合改进后的序列A和序列B再与传统PSA方法进行对比会更具有意义。
尽管如此,新方法仍然在开发与实践先进核安全分析上迈出重要一步。其整合了两种传统分析方法并且考虑了多项参数的不确定性,能够有效合理的减小堆芯熔毁概率。更重要的是,我们可以预见当面对电厂设计变更的情况时,该方法有能力敏感的量化出其对电厂造成的风险。
小幅的设计变更的风险响应变化值用传统的PSA技术是没有办法合理量化的。本节将运用CRA技术对电厂小幅设计变更的风险响应进行量化评估,并与上节中得到的额定功率下的CDF进行对比。
将堆芯功率提升5%,其他初始条件不变。基于CRA方法的PSA的改进模型,重新确定重要参数的抽样范围和分布,对电力恢复时间的区间划分也将重新定义,系统程序对于抽样算例的模拟计算也将重新进行,进而计算出功率提升下的堆芯熔毁概率,量化出功率提升下的风险变化值。
在升功率的额定工况下,由系统程序计算出堆芯温度响应以及为了避免堆芯熔毁的最晚安注时间的确定为13.968 h,该时间比功率提升前提前了16 min左右。当轴封失效发生在13.968 h以后,抽样算例属于丧失热阱为主因的序列。此时BBC和WBC的堆芯包壳温度响应曲线如图11所示,功率提升下的tmin和tmax为12.835和13.100。此时的随机性参数的抽样范围如表5所示,抽样函数不发生变化。
图11 最佳工况和最差工况的包壳峰值温度变化图(功率提升)Fig.11 Peak cladding temperature of BBC and WBC
表5 随机性参数抽样范围(功率提升)
如图11所示,不难看出当功率提升了5%时,堆芯包壳温度升高的起始时间较功率未提升时发生了提前,并且温度爬升的速度略有变快,最终导致避免堆芯熔毁的电力最晚恢复时间均有所提前。对比升功率前后的(tmin,tmax)区间,可见抽样区间变窄并且发生前移。可以预见,抽样区间的前移将会导致事故序列一定成功的概率变低,一定失败的概率变高,条件成功的概率还需要进一步确定。但无论如何,经过这些改变,功率提升的电厂风险响应得以体现,并且量化。
按照上述的参数抽样范围和抽样函数,对8个电厂状态参数和3个随机性参数进行抽样。由此,丧失热阱为主因序列的50组算例输入参数抽样完成。运用系统程序对所有算例进行模拟和计算。
与额定功率情况类似,直接计算得到SL,BS,BF和BC分支的概率。此时,轴封失效时间大于13.968 h的概率是0.007 8。
对于电力的恢复一定不能够及时抑制堆芯的熔毁的事故序列,序列3。根据新PSA模型(见图7),序列3公式如(14)。因为tmax较原始功率时有所提前,导致p(tDGr>tmax)·p(tosr>tmax)计算值较原功率的0.142 88增加到0.153 26。所以得到一定失败的序列3发生的概率为3.227 6×10-5,较原始功率增加了21%。由此直观看到功率提升对一定发生堆芯熔毁的概率的影响。
(14)
同时,由改进的PSA模型和公式(5)可以得到,功率提升情况下序列4和序列5两个序列发生的概率和1.693 4×10-6,较原始功率下的概率有所减小,主要由于(tmin,tmax)区间变小引起。
由系统程序的建模和物理过程的模拟之后,获得50个功率提升情况下对应的PCT值。散点图中如图12所示。
图12 丧失热阱为主因序列算例的50个包壳峰值温度结果散点图(功率提升)Fig.12 50 sets of PCT results in loss of heatsink dominant case(power uprate)
根据卡方检验方法得到功率提升后的PCT数据仍然满足正态分布(见图13)。假设置信度为95%,代入公式(9)和公式(10)得到总体均值和标准差为(1 469.519,203.521)。所以PCT总体数据大于1 477.6 k的概率为pce,5为48.4%。
图13 丧失热阱为主因序列算例结果直方图(功率提升)Fig.13 50 sets of PCT results in loss of heatsink dominant case(power uprate)
根据公式(11)和公式(12)计算得到序列3和序列5的堆芯熔毁概率分别为3.835×10-8/堆年和9.857×10-10/堆年。将两者相加得到改进的序列A的堆芯熔毁概率为3.934×10-8/堆年。
与功率提升之前的结果进行比较(见表6),序列3的概率升高了21.2%。因为序列3代表的是电力恢复时间过晚(t>tmax),一定来不及抑制堆芯熔毁的概率。在功率提升之后,整个核电厂的形势较没有升功率之前更为严峻,(tmin,tmax)条件成功区间由(13.017~13.447)变化成 (12.835~13.1),tmin和tmax均变小,抽样区间左移,堆芯一定熔毁对应的tmax变小导致电力恢复一定来不及恢复的概率变高。同时我们观察到,序列5的条件失效概率由原功率下的47.1%微微增加到48.4%,这同样是功率提升对电厂风险带来影响的体现。而之所以序列5最终的堆芯熔毁概率有微幅的下降,是因为条件成功区间变得狭小,发生频率有所降低,最终导致了序列5堆芯熔毁概率有所降低。然而在功率提升情况下,综合两个序列的CDF,丧失热阱为主因的事故序列的堆芯熔毁概率为3.934×10-8。对比升功率之前的堆芯熔毁概率3.28×10-8,CRA方法通过细化成功准则量化出了小幅功率提升对于电厂的风险的影响,其风险增长了19.9%。
表6 功率提升前后丧失热阱为主因的事故各序列堆芯熔毁概率对比
本文针对小幅设计变更的风险量化问题,开发了一套以RISMC方法论为指导的安全分析方法——计算风险评估(CRA)方法。该方法整合了概率论PSA技术以及确定论BEPU技术,能够综合处理随机性不确定性参数和认知性不确定性参数。特別是,对于传统PSA方法无法合理量化的小幅功率提升的风险影响,该方法展现出明显的优势。得到的具体结论如下:
(1)本文通过对真实电厂建模、对传统PSA模型关键事故进行细分,对电厂状态参数和随机性参数同时随机抽样以及对不确定性进行量化等步骤,最终计算得到以RISMC为指导、CRA方法实施下的堆芯熔毁概率。相对于传统PSA的计算结果,CRA计算方法可将堆芯熔毁概率计算值有效减小,从而挖掘出了更大的安全裕度空间。
(2)运用CRA方法量化功率小幅提升对典型压水堆全厂断电事故概率显著序列的堆芯熔毁风险的影响。在CRA的方法中,当电厂出现小幅度的功率提升时(5%),tmin和tmax会因此发生左移而增加一定失败的概率。同时,抽样范围内的条件失效概率亦随之增加,最终造成成功准则中的失败概率变高,因而功率小幅提升的风险得以合理量化(见表6)。
综上,基于RISMC方法论开发的CRA方法是一个整合两种现有安全分析技术并且综合考虑多类不确定性因素以进行风险指引的安全裕度特性分析与评估的有效途径,同时该方法能够合理量化小幅设计变更对电厂造成的风险影响。