一种压水堆大破口失水事故始发应急工况预测方法与流程

本发明属于核辐射安全技术领域，具体涉及一种压水堆大破口失水事故始发应急工况预测方法。

背景技术：

福岛事故再一次为人类敲响了警钟，再一次表明：虽然反应堆严重事故发生概率极低，但是它仍然可能发生。iaea提出，应直接根据堆芯损伤评价的结果指导防护行动，而无需依据源项估算结果。在应急情况下依据几个关键的反应堆工况参数实时、快速地预测未来的反应堆工况与事故进程，可以为应急决策提供支持，为场内和场外防护行动的实施赢得更多的时间。

由于神经网络的方法具有计算速度快这一显著优点，所以国内外学者使用该方法进行了一些事故工况预测的研究。

韩国学者sangjunpark，jinkyunpark等开发了二回路瞬态的监测系统，用以在事故早期阶段探查到瞬态的发生，确定瞬态的场景类型，并告知操作员距离汽轮机停机的剩余时间。在1000mw压水堆模拟机上的瞬态示例的数据库中，建立了所有的瞬态场景。瞬态示例展示了在没有操作员干预的情况下汽轮机停机之前核电站的行为。与此同时，从核电站的计算机上周期性地获取运行数据，并将之与瞬态示例数据库进行对比，找到与之高度相似的瞬态示例，从而确定瞬态的类型和预计的距离汽轮机停机剩余时间。对运行数据的预处理包括标准化、降低维数、确定权重因子。数据库包含上百个变量，所以很难迅速比较各个示例、及时给出结论，故使用基本组件分析方法降低变量的维度。本方法集中在预测汽轮机停机时间的计算和瞬态类型的确定上，然而本方法也可以用于在核电站中探查系统失效。主要流程为：建立二回路系统瞬态实例的数据库，与实时数据进行相似性分析，最终给出瞬态类型和距离汽轮机停机的剩余时间。详见图1。

韩国学者seunggeunkim等使用支持向量计算器，利用核电站状态变量输入在短时间内的变化，预测严重事故中可能发生的事件的发生时间。因为核电站的尺寸和复杂性一直在增加，操作员在理解核电站系统问题和缓解问题中面临着越来越大的挑战。在瞬态发生时，操作员要检查大量的仪表数据以确定异常事件的类型。由于数据量大、诊断时间短等复杂因素，操作员可能会做出错误的操作，从而加剧事故的严重性。如果在核电站发生瞬态事件，操作员会根据推荐操作规程估计异常事件的类型、尝试缓解异常事件的后果。对于操作员而言，严重事故相关事件的发生时间是非常重要的。因此，在严重事故发生时，能够提供此类信息的支持系统是大有作为的。首先，确定破口位置和破口尺寸。使用maap程序获得演练和测试数据。结果表明所开发的算法能够在loca事故中正确地确定破口位置，准确确定破口尺寸，也能够在合理的误差范围内预测主要事件的发生时间。根据这些结果，该文献提出的算法可以用于真实的核电站，因为该算法仅使用紧急停堆后短时间内的数据。如今有许多计算机代码可以进行事故相关的安全分析，但是它们的计算时间偏长，不适合作为实时支持系统。有研究使用人工智能技术处理监测数据，成功预测了事故相关事件的发生时间，包括堆芯裸露时间、堆芯温度超过1200℃的时间、压力容器失效时间。然而该研究只分析没有认为干预的情况，这种情况不符合实际。该研究将loca事故作为始发事件，分析在loca事故下的8种不同的事故序列，分析的堆型为apr1400。为了区分破口发生在热管段还是冷管段，支持向量分级(svc)使用13种核电站状态参数：破口段蒸汽发生器(sg)的压力、水位、温度；未破口段sg的压力、水位、温度；稳压器的压力和水位；堆芯水位和水温；污水坑的水位；ctmt的压力和温度。确定破口尺寸也使用上述13种电站状态参数。当衰变热大于一回路向二回路的传热量时，堆芯最大温度才可能超过1200℃，所以选择与热损失相关的参数。压力容器失效可以看作是堆芯温度上升的延伸，所以选择的参数是一样的。预测堆芯最大温度超过1200℃使用的参数有：未破损段s/g的温度，安全壳压力、温度，压力容器水位；压力容器失效时间使用的参数：未破损段s/g的温度，安全壳压力、温度，压力容器水位；安全壳失效时间使用的参数(压力超过4个大气压时，认为安全壳失效)：稳压器压力，安全壳压力和温度，压力容器水位。得出的结果主要为不同破口尺寸与主要事件发生时间的对应关系曲线；未来可以继续开展的工作主要是降低预测的误差。

神经网络方法已经广泛地应用于核工程领域中，有许多学者使用神经网络的方法进行反应堆事故进程预测的研究。韩国学者soonhopark等使用模糊神经网络预测严重事故中的压力容器水位，因为在严重事故中压力容器的水位无法测量到。该研究使用培训数据开发模型，并使用与之相互独立的测试数据进行验证；所用数据由maap4程序模拟oap1000获得。

近来，核电站的严重事故已经成为全世界关心的问题。在事故初期，核电站的主要安全参数急剧变化，使得操作员没有充足的时间做出响应。事故期间，氢气从一回路泄漏至安全壳内，氢气浓度达到一定水平后会有发生氢爆的危险，严重威胁安全壳的完整性。氢气来源于锆水反应、安全壳内的蒸汽受到辐照活化和安全壳喷淋造成内部材料的腐蚀。为了防止发生氢气爆炸，安全壳内的氢气浓度应低于4％。韩国学者dongyeongkim等使用人工智能的方法预测安全壳内氢气浓度的变化。该研究将破口位置分为热管段、冷管段以及蒸汽发生器传热管三种。破口尺寸分为210个步长；对于loca事故，破口尺寸从双端断裂尺寸的1/10000～1；对于蒸汽发生器，传热管破裂数量从1根变化到210根。破口尺寸的误差为0.4％。预测安全壳氢气浓度所用的两个输入为破口尺寸、停堆时长；使用maap4的模拟数据进行模型的开发和验证。

韩国学者mangyunna，sunhoshin等使用概率神经网络方法预测何时堆芯发生裸露、cet温度何时超过648.9℃(一般认为此时开始严重事故管理)、压力容器何时失效。热管段loca为始发事件时，预测堆芯何时裸露使用的输入有3个：安全壳温度、稳压器水位、完好段s/g温度。预测cet温度超过648.9℃所使用的输入有3个：污水池水位、压力容器水位、破口段s/g水位。模糊约束的数目为4个。预测压力容器失效的输入有3个：安全壳温度、稳压器水位、稳压器压力。冷管段loca为始发事件时，预测堆芯何时裸露使用的输入有3个：压力容器水位、污水池水位、堆芯水位。预测cet温度超过648.9℃所使用的输入有4个：压力容器水位、污水池水位、完好段s/g水位、破口段s/g水位。预测压力容器失效的输入有3个：安全壳压力、稳压器压力、破口段s/g水位。

清华大学核能与新能源技术研究院针对位于山东荣成的高温气冷堆示范电站，开发了一套核应急响应支持系统。该系统有两个任务：一是基于测量仪器的异常进行故障诊断；二是基于测量仪器的读数和操作员的行动预测事故进程。为了更好地执行这两项任务，开发了动态贝叶斯网络，在该网络中使用贝叶斯机器学习(lbp)算法进行事故诊断和预测。

综上所述，反应堆严重事故预测研究汇总，可总结为表1。

表1

针对压水堆大破口事故始发应急，现有技术中主要采用的方法是，以流体的质量守恒、和能量守恒方程为模型的基础，守恒方程的形式取决于所采用的两相流模型。但判断为气相或液相后，再采用质量守恒、和能量守恒方程计算工况评价所需参数的过程较为繁琐，需要的时间较长。

技术实现要素：

针对现有技术中存在的缺陷，本发明的目的是提供一种压水堆大破口失水事故始发应急工况预测方法，将复杂的一回路进行合理简化，对一回路进行瞬态分析，通过将事故划分为三个阶段，根据不同的阶段特点使用不同的经典公式近似估算应急工况进程，使得计算速度将大大加快，提升评价效率。

为达到以上目的，本发明采用的技术方案是：一种压水堆大破口失水事故始发应急工况预测方法，所述应用包括以下步骤：根据大破口失水事故特点将事故分为喷放、再灌水/再淹没和长期冷却三个阶段；

计算包壳温度，包壳温度变化值＝(堆芯产热量+锆水反应放热量-堆芯换热量)/燃料棒比热；

根据堆芯温度评价事故应急工况。

进一步的，堆芯换热量的计算公式为：

qtransfer＝haδt

式中，h-换热系数，w/(m²·k)；a-燃料棒外表面积，m²；δt-燃料棒温度与冷却剂温度差，k。

进一步的，在喷放和再灌水/再淹没阶段，当堆芯冷却剂与燃料棒的对流方式属于强迫对流时，所述换热系数h的计算方法包括：

计算层流强迫对流换热系数包括以下公式：

nu＝c(n)gdev

式中，c(n)-对于棒束和单管均取4.36；(z-z0)-距离流动入口的长度，1000m；dh-水力学直径；re-雷诺数；pr-普朗特数；gdev-充分发展的nusselt数；f(z)-无量纲入口长度；

计算紊流强迫对流换热系数公式为：

nu＝0.023re^0.8pr^0.4

式中，re-雷诺数；pr-普朗特数；

比较层流强迫对流换热系数和紊流强迫对流换热系数的大小，取两者最大值作为强迫对流换热的对流换热系数h。

进一步的，在长期冷却阶段，堆芯冷却剂与燃料棒的对流方式属自然对流时，所述换热系数h的计算方法公式包括以下两种，取两者最大值作为自然对流换热的对流换热系数：

nu＝0.18ra^1/4(l/dh)^-1/9

nu＝0.065ra^1/3(l/dh)^-1/9

式中：l-流道长度；ra-瑞利数；dh-表征流动的特征长度。

进一步的，所述锆水反应放热量计算方式为：

qox(t)＝qox(t0)+hrp(t)-hrp(t0)

hrp(t)＝hr(t)-hp(t)

式中：t0-参考温度，298.15k；qox(t)-温度t下锆水反应放热量，j；hr(t)-温度t下反应物焓值之和，j；hp(t)-温度t下反应物焓值之和，j；qox(t0)-锆水反应对应值为5.797×10⁶j/kgzr。

进一步的，所述堆芯产热量计算公式为：

pdecay＝n(τ)+ns1(τ)+ns2(τ)

式中：n(τ)-τ时刻剩余中子引起的裂变功率，w；ns1(τ)-τ时刻裂变产物的衰变功率，w；ns2(τ)-τ时刻中子俘获产物的衰变功率，w；

n(τ)/n(0)＝0.15exp(-0.1τ)

式中：τ-停堆时间，s；n(0)-停堆时反应堆功率，w；

式中：ns1(τ)-停堆τ秒后的裂变产物的衰变功率；a和a为系数；τ0-连续运行的时间；

式中：c和α为常数。

进一步的，当10^-1s≤τ＜10s，a＝12.05，a＝0.0639；

当10s≤τ≤1.5×10²s，a＝15.31，a＝0.1807；

当1.5×10²s≤τ＜4×10⁶s，a＝26.02，a＝0.2834；

当4×10⁶s≤τ＜2×10⁸s，a＝53.18，a＝0.3350。

进一步的，当为喷放阶段时，冷却剂温度的计算中需要的参数为堆芯压力、泄漏流量、堆芯水位；

当为再灌水/再淹没和长期冷却时，冷却剂温度的计算中需要的参数为中压安注、堆芯压力、堆芯水位。

本发明的效果在于，将复杂的一回路进行合理简化，对一回路进行瞬态分析，通过将事故划分为三个阶段，根据不同的阶段特点使用不同的经典公式近似估算应急工况进程，使得计算速度将大大加快，提升评价效率。

附图说明

图1为现有技术中诊断与预测流程示意图；

图2为压水堆大破口失水事故阶段划分原理示意图；

图3为堆芯活性区节点划分示意图；

图4为应急工况预测程序的计算流程示意图。

具体实施方式

为使本发明解决的技术问题、采用的技术方案和达到的技术效果更加清楚，下面将结合附图对本发明实施例的技术方案作进一步的详细描述。显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，均属于本发明保护的范围。

本发明提供一种压水堆大破口失水事故始发应急工况预测方法，所述方法包括以下步骤：

根据大破口失水事故特点将事故分为喷放、再灌水/再淹没和长期冷却三个阶段。参阅图2，各阶段的特点如下：

喷放阶段的特点包括：冷却剂快速泄露、一回路迅速泄压、堆芯泡份额增加以及传热恶化，堆芯温度上升；

再灌水/在淹没阶段的特点包括：中压安注开始注水、堆芯水位上升并完全淹没堆芯以及传热有所改善，堆芯温度暂时下降；

长期冷却阶段的特点包括：中压安注箱排空、衰变热加热堆芯剩余冷却剂、冷却剂汽化，从破口逸出、水位下降，堆芯逐渐裸露以及堆芯温度上升，直至熔化。

还需要说明的是，参阅图3，在本方法中将堆芯活性区沿径向划分为4个径向环，沿轴向划分为10个轴向层，共40个节点。在计算堆芯温度时，认为每个节点的温度是相同的。

还需要说明的是，参阅图4，根据堆芯压力是否大于安全壳压力来判断是否已由喷放阶段转换为再灌水/在淹没阶段；根据堆芯是否满水位判断是否已由再灌水/在淹没阶段转换为长期冷却阶段。

计算包壳温度，包壳温度变化值＝(堆芯产热量+锆水反应放热量-堆芯换热量)/燃料棒比热。

(1)堆芯换热量计算

堆芯换热量的计算公式为：

qtransfer＝haδt

式中，h-换热系数，w/(m²·k)；a-燃料棒外表面积，m²；δt-燃料棒温度与冷却剂温度差，k。

其中，在喷放和再灌水/再淹没阶段，当堆芯冷却剂与燃料棒的对流方式属于强迫对流时，所述换热系数h的计算方法包括：

计算层流强迫对流换热系数以下公式：

nu＝c(n)gdev

式中，c(n)-对于棒束和单管均取4.36；(z-z0)-距离流动入口的长度，1000m；dh-水力学直径；re-雷诺数；pr-普朗特数；gdev-充分发展的nusselt数；f(z)-无量纲入口长度；

计算紊流强迫对流换热系数公式为：

nu＝0.023re^0.8pr^0.4

式中，re-雷诺数；pr-普朗特数；

比较层流强迫对流换热系数和紊流强迫对流换热系数的大小，取两者最大值作为强迫对流换热的对流换热系数h。

在长期冷却阶段，堆芯冷却剂与燃料棒的对流方式属自然对流时，所述换热系数h的计算方法包括以下两种，取两者最大值作为自然对流换热的对流换热系数：

nu＝0.18ra^1/4(l/dh)^-1/9

nu＝0.065ra^t/3(l/dh)^-1/9

式中：l-流道长度；ra-瑞利数；dh-表征流动的特征长度。

(2)锆水反应放热量计算

若在上述阶段中发生过锆水反应，则计算锆水反应放热量。锆水反应化学方程式为：

zr+2h2o→zro2+2h2+qox

式中，qox表示反应放出的热量，j。

则，锆水反应放热量计算方式为：

qox(t)＝qox(t0)+hrp(t)-hrp(t0)

hrp(t)＝hr(t)-hp(t)

式中：t0-参考温度，298.15k；qox(t)-温度t下锆水反应放热量，j；hr(t)-温度t下反应物焓值之和，j；hp(t)-温度t下反应物焓值之和，j；qox(t0)-锆水反应对应值为5.797×10⁶j/kgzr。

(3)堆芯产热量计算

堆芯产热量计算公式为：

pdecay＝n(τ)+ns1(τ)+ns2(τ)

式中：n(τ)-τ时刻剩余中子引起的裂变功率，w；ns1(τ)-τ时刻裂变产物的衰变功率，w；ns2(τ)-τ时刻中子俘获产物的衰变功率，w；

n(τ)/n(0)＝0.15exp(-0.1τ)

式中：τ-停堆时间，s；n(0)-停堆时反应堆功率，w；

式中：ns1(τ)-停堆τ秒后的裂变产物的衰变功率；a和a为系数；τ0-连续运行的时间。根据不同的堆停时间，a和a取不同的值，具体查下表2可得。

表2

中子俘获产物的衰变功率ns2(τ)可由下式计算：

对于低富集度铀作燃料的压水堆，可取c＝0.6，α＝0.2。在其他实施例中，中子俘获产物的衰变功率ns2(τ)在上式所得数值基础上再乘以系数1.1。

根据堆芯温度评价事故应急工况。

根据包壳温度可判断堆芯工况，判断依据见表3。原始堆芯几何形状可在数分钟至数小时的时间内逐渐丧失，对于轻水堆，温度范围在1500k～3000k。对于低压事故序列，堆芯温度达到1000k～1200k时，锆包壳开始肿胀、破裂；对于高压事故序列，包壳在低温下坍塌到燃料上，锆合金包壳失效时间可能推迟，堆芯温度达到1500k时才会开始肿胀、破裂。

就严重事故而言，因堆芯材料氧化而放出的热量可能超出衰变热。金属锆氧化释热率在1500k时与衰变热相当，在1800k以上的温度下可能达到衰变热的10倍。

堆芯材料的整体熔化取决于每种堆芯材料的熔化温度，但由于堆芯内不同材料在相互紧密接触的情况下发生化学反应，其整体熔化会受到所形成低熔点合金的影响。由于因科镍合金、不锈钢和包壳的剧烈相互作用，定位格架将因低熔点合金的形成二在相对低的温度下而受损。对于典型压水堆设计，在1500k～1700k温度范围内，fe-zr、ag-zr等之间的化学反应会导致堆芯构件的早期液化与迁移。

锆合金的熔点一般介于2000k～2250k之间，包壳温度一旦超过2000k，锆合金包壳就会熔化，而且在某些情况下会流动到堆芯下部区域。熔融的锆会部分溶解uo2芯块和包壳外的zro2，使其在远低于自身熔点的温度下液化。

表3

在堆芯换热量的计算中运用到了δt(燃料棒温度与冷却剂温度差)。其中关于冷却剂温度的计算，根据不同的阶段特点，参阅图4计算所需的参数也不相同。

当为喷放阶段时，冷却剂温度的计算中需要的参数为堆芯压力、泄漏流量、堆芯水位；当为再灌水/再淹没和长期冷却时，冷却剂温度的计算中需要的参数为中压安注、堆芯压力、堆芯水位。

(1)中压安注

利用伯努利方程计算中压安注流量，当中压安注箱中的水耗尽或中压安注箱的压力小于等于一回路压力时，流量为0。

(2)堆芯压力

堆芯压力为冷却剂温度对应的饱和压力，使用下式计算：

式中：ps-水蒸汽的饱和压力，kpa；t-水蒸汽温度，k；a、b、c为常数，查下表4可得。

表4

(3)泄漏流量计算

通过流入和流出一回路的水质量平衡来计算泄漏流量。流入流出一回路的各流量、一回路压力温度等，使用核电厂的实时测量数据。

泄漏流量使用下式计算：

式中，wb-泄漏流量，kg/s；win-冷却剂流入一回路的流量，kg/s；wout-冷却剂流出一回路的流量kg/s；-计算时间间隔内，一回路冷却剂因温度压力变化而产生的质量变化，kg/s。

冷却剂流入一回路的流量使用下式计算：

win＝wsi+wfeed+wseal+wspray

式中，wsi-安注流量，kg/s；wfeed-化容补水流量，kg/s；wseal-主泵轴封注入流量，kg/s；wspray-稳压器辅助喷淋流量，kg/s。

冷却剂流出一回路的流量使用下式计算：

wout＝wletdown+wsl

式中，wletdown-一回路下泄流量，kg/s；wsl-主泵轴封回水流量，kg/s。

式中，-一回路冷却剂除稳压器外的平均密度，kg/m³；v-一回路冷却剂系统除稳压器外的体积，m³；ρl_preρl_pre-稳压器液态冷却剂密度，kg/m³；s-稳压器横截面积，m²；hpre-稳压器液位，m；ρg_pre-稳压器气态冷却剂密度，kg/m³；hpre-稳压器总高度，m。

(4)堆芯水位

初始堆芯冷却剂总质量减去冷却剂泄漏总质量即为堆芯冷却剂剩余质量(与泄漏量相比蒸发量很小可以忽略)，剩余冷却剂除以冷却剂密度得到剩余冷却剂体积，剩余冷却剂体积除以压力容器横截面积即可得到堆芯水位，即：

式中：mtotal-初始堆芯冷却剂总质量，kg；mleak-冷却剂泄漏总质量，kg；ρcoolant-冷却剂密度，kg/m³；srpv-压力容器横截面积，m²。

(5)冷却剂温度计算

气态冷却剂与裸露的燃料元件换热，液态冷却剂与未裸露的燃料元件换热，温度上升均取决于换热量及其比热。由冷却剂温度初始值、每时间步长的温度变化，可知未来冷却剂随时间变化。

δtcoolant_g＝qtran_g/cp_g

式中：δtcoolant_g-气态冷却剂温度变化；qtran_g-气态冷却剂与堆芯换热量；cp_g-气态冷却剂比热。

δtcoolant_l＝qtran_l/cp_l

式中：δtcoolant_l-液态冷却剂温度变化；qtran_l-液态冷却剂与堆芯换热量；cp_l-液态冷却剂比热。

区别于现有技术，本发明提供的一种压水堆大破口失水事故始发应急工况预测方法，将复杂的一回路进行合理简化，对一回路进行瞬态分析，通过将事故划分为三个阶段，根据不同的阶段特点使用不同的经典公式近似估算应急工况进程，使得计算速度将大大加快，提升评价效率。

本领域技术人员应该明白，本发明所述的方法并不限于具体实施方式中所述的实施例，上面的具体描述只是为了解释本发明的目的，并非用于限制本发明。本领域技术人员根据本发明的技术方案得出其他的实施方式，同样属于本发明的技术创新范围，本发明的保护范围由权利要求及其等同物限定。