示范快堆一回路主管道断裂事故最佳估算分析研究

叶尚尚，刘一哲，杨红义，杨 军，王晓坤，齐少璞

(中国原子能科学研究院，北京 102413)

与安全裕量有关的研究一直是反应堆安全设计与安全分析中的重点和难点问题。核电厂的安全评价方法主要有两种：一种是依据设计基准事故的确定论评价法；另一种是概率安全评价(PSA)法。确定论方法是核电厂发展史上长期使用的方法，广泛应用于世界各国的核电站安全分析中，其基本思想是依据反应堆纵深防御的原则，除在设计上尽可能安全可靠外，还设置了多重专设安全设施，以便在发生最大假想事故情况下，依靠安全设施，将事故后果减轻至最轻程度。在确定安全设施的种类、容量和响应速度时，需要一个参考的假想事故作为设计基础，并将这一事故看作包络性事故，认为所设置的安全设施若能防范这一事故，就必能防范其他各种事故。PSA方法是应用概率风险理论对核电厂安全性进行评价，认为核电厂事故是随机事件，引起核电厂事故的潜在因素很多，核电厂的安全性应由全部潜在事故的数学期望值表示。

随着国际核电行业的迅猛发展和核安全要求的不断更新，国际上确定了4种用于执照申请的安全分析方法，分别为保守性方法、组合方法、最佳估算分析方法和风险指引方法。其中，最佳估算分析方法能提供更接近反应堆行为的信息、区分大部分相关安全问题以及较准确地获得安全裕量，是目前现实可行的分析方法，逐渐成为事故安全评审的主要发展趋势。

示范快堆CFR600是一座池式钠冷快堆，采用钠-钠-水三回路两环路布置，一回路每环路设置1台主循环泵，每台泵出口至大栅板联箱有2条压力管线。CFR600在发生一回路主管道断裂、主泵卡轴等一回路失流事故时，堆芯流量迅速降低，包壳温度出现峰值，在叠加极端保守假设的情况下，会造成堆芯损伤。为更确切地确定事故分析安全裕量，并量化不确定性结果，降低采用保守方法进行事故分析对电厂经济性与灵活性的限制，本文采用最佳估算加不确定性(BEPU)分析方法对CFR600一回路主管道断裂事故进行分析。

1 最佳估算分析方法

最佳估算是指尽可能真实地反映物理现象，不人为引入悲观性，使用最佳估算程序并对结果进行不确定性分析[1]。最佳估算的核心技术主要包括最佳估算程序和不确定性分析技术。

1.1 最佳估算模型和程序

最佳估算程序的核心是反映重要物理现象的计算模型，因此应分析池式钠冷快堆主热传输系统热工水力重要现象，基于重要现象开发能真实反映系统物理瞬态特性的热工水力模型，该模型能满足所需分析瞬态的要求，并充分模拟其中的重要物理现象。开发的模型需要采用试验数据或基准例题对模型的适宜性进行评估。

池式钠冷快堆主热传输一回路系统采用一体化布置，具有典型的特殊现象，包括钠池热分层现象、中间热交换器一次侧倒流现象、主管道断裂喷放现象以及一回路主泵惰转现象等。主热传输二回路系统采用回路式布置，系统中的重要现象包括单相/两相界面滑移现象、水/汽两相换热现象(指蒸汽发生器传热管内介质)以及二回路自然循环流动现象等。基于主热传输系统重要现象分析，分别建立主热传输系统中关键部件的最佳估算模型。

1)基于时间常数的主泵惰转模型[2]

当核电厂发生失电事故时，由于飞轮和管路内冷却剂流动的惯性，主循环泵仍以瞬变转速持续转动。该瞬态过程主要分为2个阶段：第1阶段，在瞬变开始时，主循环泵的惯性压头较重力压头大得多，前者与泵的转动惯量有关，后者与主泵所在回路的流动惯性有关；第2阶段，在惰转后期，主泵的转速逐渐下降为0，其惯性压头逐渐消失，流体完全依靠流动惯性驱动，即稳态自然循环。

传统模型中，主泵惰转高转速阶段和低转速阶段分别考虑了泵的惯性压头(即转速的2次方项占主导地位)和回路阻力特性(转速较低时，一般认为克服冷却剂流动所需的水力学力矩与转子的摩擦力矩之和为常数)，事实上，在泵的整个惰转过程中均需考虑上述2项。基于泵的转矩方程，经过合理分析简化，推导出基于时间常数的泵惰转特性曲线计算模型，无需明确泵的转动惯量和回路系统的阻力特性，只需输入2个时间常数即可，如半转速时间和到零时间，并采用主循环泵的试验数据进行验证。

2)基于射流/羽流理论的热分层模型[3]

建立了基于羽流/射流理论的钠池热分层效应模型，以解决系统程序零维或一维模型较难分析出钠池热分层现象对系统热工水力特性的影响。

从喷嘴喷射出的流体的运动过程如图1所示，随着流体速度的降低，流体高度减小，横截面积增大。

图1 射流发展示意图Fig.1 Schematic diagram of jet development

射流能达到的高度与流体射出的初始速度、初始密度和环境介质密度相关，推荐的用于计算射流高度的经验关系式如下：

(1)

式中：Zjet为射流高度，m；Fr0为无量纲数；r0为射流入口半径，m。

Fr0的计算式为：

(2)

式中：Vin为射流入口速度，m/s；ρin、ρB分别为入口密度和环境介质密度，kg/m3。

根据射流理论的思想，将热池空间划分成上下2个区域，区域内部认为完全搅浑，区域分界面的高度一部分由射流高度决定，同时，考虑当射流高度低于出口管嘴时由进入流体完成的填充高度带来的影响。基于射流/羽流理论的钠池模型示意图如图2所示。

图2 基于射流/羽流理论的钠池模型示意图Fig.2 Schematic diagram of sodium pool model based on jet/plume theory

区域分界面高度Z(t)的计算式如下：

Z(t)=Zjet(t)Zjet(t)≥Hext

(3)

Zjet(t)

(4)

式中，Zf(t)为填充高度，m，其计算公式为：

(5)

当区域A、B的界面确定后，其能量守恒方程为：

hNa-mAmA(TmA-TPoolA)+xρPoolBANa·

|dZ(t)|cpPoolB(TPoolB-TPoolA)

(6)

hNa-NaANa(TPoolA-TPoolB)+hNa-mAmB·

(TmB-TPoolB)+(1-x)ρPoolAANa·

|dZ(t)|cpPoolA(TPoolA-TPoolB)

(7)

式中：M为质量，kg；T为温度，℃；t为时间，s；A为结构件与流体的换热面积，m2；Win为入口质量流量，kg/s；hNa-Na和hNa-m分别为区域A、B之间的换热系数及区域与金属结构件之间的换热系数，J/(m2·℃)；下标m和Pool分别指结构件和流体，PoolA、PoolB分别代表区域A和区域B；x表示因界面高度变化对区域A、B质量造成的影响，其特征量为：

(8)

与区域A、B相对应的金属结构件的能量方程为：

(9)

(10)

出口温度根据分界面的高度可表达成区域A、B温度的函数：

Tout=yTPoolA+(1-y)TPoolB

(11)

(12)

3)一回路主管道断裂模型[4]

池式钠冷快堆一回路采用池式布置，因此流网较复杂，简化后的一回路系统流网示意图如图3所示。其中，方框内的编号代表节点，连线为等效的管道。

一回路系统的流网采用特征线方法进行求解。采用图3所示管网进行一回路主管道断裂事故模拟，在发生断裂的压力管上设置一跳脱阀，同时在栅板联箱和冷池以及泵出口和冷池之间添加管道，管道的总长度和直径与压力管相当，管道上布置破口阀。正常运行时，跳脱阀打开，破口阀1、2处于关闭状态，当发生一回路主管道断裂事故时，跳脱阀瞬间关闭并阻止该管道中的冷却剂流动，破口阀1、2瞬间打开，以模拟破裂口的流量喷放现象。

图3 一回路系统流网示意图Fig.3 Flow network diagram of primary loop system

跳脱阀和破口阀完全打开时的流通面积与主管道(单根压力管)相同，完全关闭时接近0。跳脱阀完全打开时的阻力与对应长度的压力管阻力相同。对于完全打开后的破口阀，由于压力管流通面积相对于钠池非常小，破口阻力简化成喷向无限大空间进行计算。跳脱阀和破口阀完全关闭后的阻力均近似无限大。

由于主管道断裂现象在实际反应堆中尚未发生，且由于试验具有破坏性，因此目前尚无数据库支持开展模型验证。

4)中间热交换器热工水力模型[5]

钠冷快堆的中间热交换器(IHX)壳侧流体为一回路高温冷却剂，二次侧为二回路载热剂，两侧流体逆向流动换热。通过开发IHX热工水力模型，实现对IHX一次侧产生的倒流现象进行模拟。IHX两侧流体和传热管壁温分布计算控制方程如式(13)～(15)所示，网格划分示意图如图4所示。

图4 中间热交换器节点划分示意图Fig.4 Schematic diagram of node division of IHX

一次侧流体能量方程：

(13)

二次侧流体能量方程：

(14)

管壁传热方程：

(15)

式中：下标p表示壳程(一次侧)流体；s表示管程(二次侧)流体；m表示管壁。

5)直流式蒸汽发生器滑移网格模型[6]

钠冷快堆采用一次通过式直流蒸汽发生器(OTSG)将给水从过冷态加热至两相，并最终获取高温高压的过热蒸汽。瞬态过程中，随着边界条件的变化，相界面发生移动。为描述瞬态过程中蒸汽发生器的热工特性，避免固定网格模型容易出现的换热分区边界点阶跃和计算不稳定的缺点，采用传热区边界在瞬态过程中可移动的滑移网格模型(图5)，可有效解决瞬态过程中固定网格模型存在的界面跳跃现象。

图5 蒸汽发生器滑移网格模型简化示意图Fig.5 Simplified schematic diagram of steam generator slip grid model

基于Leibnitz理论将包含时间项的微分方程表示如下：

(16)

式中：f可为ρ、ρH，ρ为密度，kg/m3，H为焓，J/(Kg·℃)；fi=f(Zi，t)；Zi为第i控制体与给水入口的距离，m。

式(16)等于号右边第1项可表达为：

(17)

(18)

据此可得到两侧流体及管壁中含时间项方程的积分形式：

(19)

水/蒸汽两相换热系数以及传热恶化点的计算参见文献[7-8]。

基于FR-Sdaso程序[9]，通过嵌入上述关键核心模型，构成最佳估算分析程序。文献[3]采用实验数据、核动力厂运行数据、国际基准例题、成熟软件数据等对部件级核心模型和系统整体模型进行了充分的验证与确认，模型或其组成子块具有合理预测实验现象的能力以及获得与性能指标相关的预期结果的固有能力，可用于分析相似的瞬态工况[3]。

1.2 不确定性的定量化计算

目前，工程应用中评估不确定性的方法主要有7种[10]：极限值法、敏感性分析法、响应面法、快速概率积分法、基于数据库的方法、蒙特卡罗法、容忍限法。其中，容忍限法中的非参数统计可大量减少抽样工作。

非参数统计法通常不用基于某一特定的分布，而采用任意分布(或与分布无关)。当一个假设的分布被适当的实验(未知的)排除，通过随机抽样问题中的特征来决定容忍限是可能的。较早提出这种方法的是Wilks。Wilks方法[11]是一种有序的容忍限值方法，能对未知分布的随机样本建立容忍置信区间。其数学思想可描述为：对于任意变量x，建立1个容忍置信区间(L，U)，使得x所有取值至少有γ份额落在此区间的概率(置信度)为β，即：

(20)

式中：f(x)为变量x的概率密度函数，大多情况下是未知的；P为概率。

从安全分析的角度，若此置信区间上界U在所要求的限值以下，即认为是安全的。

根据Wilks的研究，对于连续的一个总体，其分布P来自于随机抽样的两个既定统计之间的总体份额，是独立于总体的。即在来源于一个随机样本中的两个有序统计量之间，整个总体中的份额是独立的，只是一个特殊的选择有序统计函数。在这种情况下，容限通过既定统计(最高、第2高等)给出。当容限通过最大既定统计给出时，需要的最小计算次数由Wilks公式给出。对于单侧和双侧容忍区间，Wilks经典公式[12-13]为：

(21)

Wilks公式仅针对1个输出量，Guba等[13]发展了Wilks理论，提出对于P个输出量(P>1)，最小计算次数由下式确定：

(22)

对于单侧容忍限，若β=γ=0.95，则N最小值为59。不同置信水平下的最小计算次数列于表1。

表1 Wilks公式确定的最小计算次数Table 1 Minimum calculation times determined by Wilks method

美国核管会(NRC)认为在95%置信度下，95%的概率已足以满足高概率的准则要求[14]。当β=95%、γ=95%时，N=59，即只需抽样59次，便可保证所有可能的结果中有95%的份额落在这59次计算结果最大值以下的概率为95%。若计算124组，则计算结果中的第3大值即为95%/95%上限。非参数抽样使得抽样次数大幅减小，计算效率得到提高，但其给出的不是目标参数的概率密度分布，而是目标参数值的置信概率和份额。

2 主管道断裂事故基准工况分析

2.1 瞬态特点及重要现象分析

CFR600一回路主管道示意图如图6所示。反应堆在额定功率下运行时，由于压力管的疲劳腐蚀、小破口、焊接缺陷或其他意外情况，造成4根压力管中的1根瞬时双端断裂且断口完全错开，从而使大量冷却剂快速从2个断裂口喷放流失，这是该事故的主要特征[15]。

图6 CFR600一回路主管道示意图Fig.6 Schematic diagram of primary loop main pipe of CFR600

事故发生后，由于反应堆主热传输回路失去平衡，造成2台一回路主泵的流量突增，并伴随着通过堆芯的冷却剂流量骤减，从而使反应堆出口温度急剧上升，反应堆会因核功率与堆芯流量之比超过整定值的保护信号而实施紧急停堆。反应堆紧急停堆信号发出后，控制棒下落，二回路主泵按照自然惰转规律惰转至停运。考虑到应急电源的作用，一回路泵按照自然惰转规律惰转至低转速。但作为单一故障，认为供给完好环路一回路泵的应急电源失效，所以完好环路一回路泵惰转至停运。此后，堆芯保持较小的流量，堆芯剩余发热由事故余热排出系统逐渐排出。

2.2 瞬态基准工况分析

在开展不确定性分析前，对该事故的基准工况进行分析，此时不考虑输入参数和边界条件的不确定性。基准工况用于分析电厂预期的真实情况，其分析结果可用于电厂设计、模型评价以及瞬态物理现象分析。

假设0时刻发生一回路1根主管道瞬时双端断裂，该事故下反应堆功率和堆芯相对流量变化示于图7，堆芯流量急剧降低，反应堆因功率流量比高，保护信号触发保护停堆。

图7 主管道断裂事故下核功率和堆芯流量变化Fig.7 Variation of nuclear power and core flow rate in main pipe break accident

当发生1根压力管双端断裂时，由于泵出口和大栅板联箱压力较冷池压力高，大量一回路冷却剂从断裂口流出，由大栅板联箱和泵出口喷放至冷池的流量示于图8，2个环路完好压力管的流量变化示于图9。

图8 主管道断裂事故下断裂口喷放流量变化Fig.8 Discharge flow rate change in main pipe break accident

图9 主管道断裂事故下完好压力管流量变化Fig.9 Flow rate change of intact pipe in main pipe break accident

1根压力管断裂后，一回路系统阻力特性发生改变，并联的2个环路流量特性发生变化，故障环路一回路流量迅速降低，而完好环路流量在短时间内升高。由于故障环路冷池压力升高，使得故障环路冷池内的液态钠打入IHX一次侧出口窗，在IHX中发生倒流，从入口窗流出进入热池。在反应堆停堆后，2台一回路主泵开始惰转，通过完好压力管的流量开始下降，一回路2个环路的流量变化示于图10。在约50 s时，由于完好环路泵产生的压头和故障环路相当，故障环路一回路流量由负变为正。在约55 s时，由于完好环路一回路泵惰转至接近零转速，从故障环路完好压力管打入大栅板联箱的钠通过完好环路IHX倒流进入热池。

图10 主管道断裂事故下一回路2个环路流量变化Fig.10 Flow rate change of two circuits of primary loop in main pipe break accident

事故下堆芯进出口温度的变化示于图11。由于冷池具有较大的钠装量，在事故发生后500 s内，冷池平均温度基本不发生变化，但由于IHX排热能力降低，导致热池热钠进入冷池，因此，堆芯入口温度在1 500 s时升至约500 ℃。堆芯出口温度在瞬态过程中出现3个峰值，第1峰值的出现是由于事故发生后堆芯流量急剧降低，功率流量比失配；完好环路一回路主循环泵的停运导致堆芯流量与功率出现第2次失配，这形成了堆芯出口温度的第2峰值；第3峰值的出现是由于堆芯入口温度的升高抬升了堆芯出口温度。事故过程中，堆芯出口钠温最高值出现在第1峰，堆芯出口平均温度最高值为668.5 ℃，组件出口钠温最高值为779.6 ℃，也在第1峰。

图11 主管道断裂事故下堆芯进出口温度变化Fig.11 Core inlet and outlet temperatures in main pipe break accident

事故下燃料包壳最高温度的变化示于图12，包壳温度的变化趋势与堆芯出口钠温的变化趋势基本一致，包壳温度最高值(806.1 ℃)出现在第1峰。

图12 主管道断裂事故下包壳最高温度变化Fig.12 Maximum temperature change of cladding in main pipe break accident

2.3 基准工况与保守计算结果的对比

事故基准工况分析不考虑电厂边界条件和初始条件的不确定性，有利于揭示事故瞬态过程中的物理现象。保守分析由于采用保守的边界条件和初始输入参数，因此得到的结果偏保守。

保守方法与基准工况下组件出口最高钠温与燃料包壳最高温度的变化示于图13。保守方法计算得到的组件出口最高钠温第1峰值达905.8 ℃，较基准工况高126.2 ℃。保守方法计算得到的包壳最高温度第1峰值达943.4 ℃，较基准工况高137.3 ℃，超出包壳破损限值(900 ℃)。

图13 主管道断裂事故下组件出口最高钠温和包壳最高温度对比Fig.13 Maximum temperature of assembly outlet and cladding of main pipe break accident

3 结果与讨论

3.1 不确定性输入参数的选取和确定

重要输入参数的概率密度分布是后续抽样的条件。这些参数的概率分布应通过试验数据或核电厂获得，具体方法[14]包括：1)从已有试验数据或程序分析文献中查找范围及其分布；2)拟合单项或整体试验中可用的试验数据，如果试验数据充足，则可用经典的统计学方法将数据转化为输入参数的范围和分布，否则可用贝叶斯方法从少量数据中推导范围和分布；3)对于核电厂运行条件，将核电厂测量值转化为输入参数的分布；4)对与热工水力模型相关的重要输入参数，确定其分布相对复杂，需用模型变量描述对应的物理过程，如用临界流的试验值与程序计算值之比作为描述程序预测临界流能力的模型变量，其分布函数由程序模拟若干临界流试验的结果分布确定。为此，有必要证明模拟具体物理过程的程序是正确的、不会引起重大偏差。

高质量试验数据以及对输入参数相互影响认识的缺乏，使得确定输入参数概率分布成为难点。对部分输入参数，若其值能严格控制，则在分析时可用其整定值，而无需考虑其不确定度分布；若参数随运行历史而变化或在事故时无法确定，则可考虑用保守值代替；若分析输入参数概率分布的代价较大，则用包络值或假设为均匀分布；若信息缺少，则通常假设各输入参数相互独立。

由于钠冷快堆运行经验较少，很多参数的分布特性尚未完全掌握，因此根据瞬态和事故分析经验，并考虑测量仪表的精度范围，选取影响池式钠冷快堆典型事故目标参数的输入参数，如表2所列。

表2 池式钠冷快堆失流事故典型不确定性输入参数Table 2 Typical uncertainty input parameter for loss of flow accident for pool SFR

续表2

3.2 不确定性计算分析

对表2中19个输入参数的124组参数进行简单随机抽样，抽样结果经归一化处理后其分布示于图14。采用最佳估算程序，以124组随机抽样参数为输入，对一回路主管道断裂事故进行瞬态计算，得到不同工况下燃料包壳最高温度的变化，如图15所示。可见，124组工况的包壳最高温度变化特性类似，由于第3峰值不会超出第1峰值，且限于计算资源的约束，图中仅示出了瞬态中前2个峰值的变化。

图14 主管道断裂事故下输入参数归一化分布Fig.14 Normalized distribution of input parameter in main pipe break accident

图15 主管道断裂事故下124组包壳最高温度Fig.15 124 maximum cladding temperatures of main pipe break accident

对计算得出的124组工况包壳最高温度(PCT)进行统计分析，得到PCT散点图和频数分布直方图，如图16所示。根据Wilks理论，124组PCT中的第3大者851.6 ℃即为95%/95%PCT的上限，即认为有95%的概率落入以上限为851.6 ℃的95%的区间内，该数值相较于包壳破损验收准则900 ℃有约48.4 ℃的裕量。

图16 主管道断裂事故下包壳最高温度统计Fig.16 Statistics of maximum cladding temperatures of main pipe break accident

最佳估算分析结果与保守分析结果的对比示于图17。不确定性计算得到的95%/95%PCT上限值介于保守方法和基准工况之间，与两者的差值分别为91.8 ℃和45.5 ℃。图17中同时给出了不确定性计算的下限，包壳最高温度为773.5 ℃，下限与上限一同构成了计算结果的不确定性范围。

图17 主管道断裂事故下包壳最高温度最佳估算与保守方法计算结果对比Fig.17 Maximum cladding temperature comparison in main pipe break accident under conservative and BEPU method

3.3 输入参数敏感性分析

敏感性分析是通过研究输入参数对目标参数的影响，识别对目标参数有重要影响的关键参数。参数敏感性分析有利于排出不太敏感或重要性相对较低的参数，从而降低参数识别的维数，提高识别的效率；或根据分析结果，进一步确定对目标参数具有重要影响的若干参数。

敏感性度量可使用Pearson、Spearman秩相关系数或Pearson、Spearman偏相关系数及标准回归系数等[16-18]。非参数抽样通常广泛使用Spearman秩相关系数ρ，其表达式为：

ρ=

(22)

式中：RXi为Xi在变量X中数值的排序；RYi为Yi在变量Y中数值的排序；n为抽样次数。

Spearman秩相关系数绝对值的大小表明各输入参数的不确定度对目标参数的敏感程度，符号表示正负相关性，即正号表示目标参数随输入参数的增大而增大；负号表示目标参数随输入参数的增大而减小。Spearman秩相关系数数值大小表征输入参数对目标参数的影响程度。

图18为一回路主管道断裂事故中PCT对输入参数的Spearman秩相关系数，实际上也是输入参数重要度。其中，功率径向分布、包壳热点因子以及主管道断裂喷放系数对PCT影响较大，在工程设计上可重点优化上述参数。

图18 主管道断裂事故下输入参数的Spearman秩相关系数Fig.18 Spearman rank correlation coefficient of input parameter for main pipe break accident

4 结论

最佳估算分析方法能较准确地获得安全裕量，逐渐成为事故安全评审的主要发展趋势。本文采用最佳估算模型和Wilks方法对CFR600一回路主管道断裂事故进行了最佳估算分析和不确定性定量计算。结果表明，一回路主管道断裂后，大量冷却剂快速从两个断裂口喷放流失，故障环路出现明显的倒流，两环路的输热特性形成明显的不对称现象。一回路主管道断裂事故下包壳最高温度95%/95%的上限为851.6 ℃，较保守分析结果有约91.8 ℃的裕量。采用Spearman秩相关系数对输入参数的敏感性进行分析，得到了影响事故后果的重要输入参数，为设计优化提供了参考和依据。