压水堆核电厂严重事故分析方法与流程

本发明属于反应堆严重事故现象计算领域，具体涉及压水堆核电厂严重事故分析方法。

背景技术：

核电厂严重事故可能带来大量的放射性物质泄漏，从而造成极大的经济损失和恶劣的社会后果，人类和平利用核能的道路上曾发生过三次比较严重的核电站事故：即美国的三哩岛事故(1979年)、前苏联的切尔诺贝利事故(1986年)和日本福岛核电站事故(2011年)。核事故将会对本世纪以来全球范围内的核电复苏抹上一层浓重的阴影。

近几年我国核电发展迅速，为了改变我国现有的能源结构，增强能源安全性、缓解对环境的压力、实现减排目标，国家制定了庞大的核电中长期发展规划。在我国核电快速发展时日本福岛核事故对我国核电安全性敲响了警钟。事故发生后，我国政府高度重视我国核电的安全问题，并明确提出了‘要对当前在役的核电厂进行安全评估’和‘要在确保安全的前提下积极发展我国核电事业’，这表明我国政府对核安全的高度重视和加强核电站严重事故研究及评估的重要性和必要性。

最近几年，国内在严重事故研究方面开始投入较多的人力和物力，并取得了很多研究成果，但是与核工业发达国家相比，我国对严重事故的研究相对不足。相对于国际上的这些较为成熟的大型商业软件，我国在严重事故程序分析开发方面缺少成熟的相关的研究。目前我国尚没有具有自主知识产权的大型先进压水堆严重事故综合分析软件，严重事故综合分析平台的建设，难度高、耗资大。而开发大型先进压水堆严重事故综合分析平台对我国核电安全分析有着重要的意义，具有自主知识产权的严重事故综合分析平台可以用于我国大型压水堆核电站严重事故分析和评审，同时开发严重事故相关的分析软件，对严重事故机理分析、事故仿真和预测、事故防范演练、以及事故预防措施研究和事故后缓解措施管理都有重要作用。

技术实现要素：

为解决上述问题，本发明在已经积累的相当的知识、经验和吸收国际上近30年严重事故最新研究成果的基础上，通过整合现有的资源和人员力量，同时基于现有严重事故分析程序中的成熟模型，开发出压水堆核电厂严重事故分析方法，分析严重事故早期堆芯行为，堆芯过热氧化熔化过程等堆内复杂的物理化学过程，堆芯形成碎片床的过程，及压力容器内熔融物行为分析。

为了到达上述目的，本发明采用如下技术方案：

压水堆核电厂严重事故分析方法，包括如下步骤：

步骤1：进行严重事故早期堆芯行为特性及燃料组件应力-应变特性计算，计算得到各堆芯节点的分布温度、堆芯包壳氧化量、氢气产生质量以及燃料包壳应力应变的变化率，早期堆芯行为特性及燃料组件应力-应变特性计算具体包含以下内容：

1)读取所有事故相态中堆芯的热工水力参数和响应信息；

2)对计算堆芯进行几何建模：根据计算需求对计算区域进行剖分，给出了一个7×13的堆芯节点划分，有三个非燃料层：两个在底部，代表堆芯支撑板和下部固定板以及下部气体腔室；一个在顶部，代表上部气体腔室和上部固定板；其他的层轴向长度相同，代表堆芯活性区；根据计算需求对计算区域进行剖分，在每一个节点内生成初始状态参数，对每个计算子区域的材料中的各状态进行赋初值；

3)基于2)所得到的初始状态参数和各材料的状态信息，利用龙格库塔法进行计算得到堆芯各节点的温度分布tn，所述的堆芯各节点温度分布tn的计算如公式(1)：

式中：

qn——流过堆芯节点的热流量；

wst——水质量流量；

hst——在冷却剂温度为tn，分压为时的焓值；

——冷却剂饱和分压；

4)基于3)所得到堆芯各节点温度分布，计算得到堆芯包壳氧化量所述的堆芯包壳氧化量的计算如公式(2)：

式中：

t——堆芯节点温度即包壳，控制棒或者水棒，燃料包壳的温度；

x——氧化物厚度；

ρzr——锆合金密度；

a——氧化系数；

b——温度修正系数；

r——理想气体常数；

5)基于4)所得到堆芯包壳氧化量，计算包壳氧化过程中的氢气产生质量，所述的氢气产生质量如反应式(5)所示，所述的氢气产生质量计算公式如(6)所示：

zr+2h2o→zro2+2h2+δhzr(5)

式中：

δhzr——每摩尔的锆合金反应产生的热量；

hst——蒸汽入口焓值；

——金属温度下的氢气焓值；

——氢气的质量流量；

wst,rct——蒸汽的质量流量；

qrct——反应产生的化学热；

——锆合金消耗的摩尔速率；

6)基于3)和5)所分别得到的堆芯各节点的温度分布和堆芯包壳氧化量，计算得到各堆芯节点的燃料包壳应力应变的变化率，所述的燃料包壳应力应变变化率的计算如公式(7)、(8)：

式中：

σr，σθ，σz——径向应力，周向应力，轴向应力；

εr，εθ，εz——径向应变，周向应变，轴向应变；

——径向塑性应变量，周向塑性应变量，轴向塑性应变量；

——径向塑性应变增量，周向塑性应变增量，轴向塑性应变增量；

——径向热变形应变量，周向热变形应变量，轴向热变形应变量；

——径向弹性变形应变量，周向弹性变形应变量，轴向弹性变形应变量；

e——杨氏模量；

v——泊松比；

αr，αθ，αz——径向热膨胀率，周向热膨胀率，轴向热膨胀率；

——膨胀导致的热应变；

boz——boltzmann常数；

k1,k2,k3,ed——经验常数；

t——温度；

由于燃料棒的长度远大于其半径尺寸，计算时将其视为一维平面应力问题进行处理；

7)利用6)计算得到各堆芯节点的燃料包壳应力应变的变化率进行下一步早期堆芯行为特性及燃料组件应力-应变特性计算，即重复3)、4)、5)、6)的计算过程，直到计算结束；

步骤2：基于步骤1中计算得到的各堆芯节点下的温度分布、堆芯包壳氧化量、氢气产生质量以及燃料包壳应力应变的变化率进行堆芯熔化计算，堆芯熔化计算相比于早期堆芯行为特性及燃料组件应力-应变特性计算在时间上是一个顺向的过程，堆芯熔化计算具体包含以下内容：

1)读取早期堆芯行为特性及燃料组件应力-应变特性计算得到的几何区域信息以及堆芯节点各状态参数；

2)利用各堆芯节点下的温度分布、堆芯包壳氧化量、氢气产生质量以及燃料包壳应力应变的变化率，对燃料组件各参数进行初始化计算，首先计算各燃料组件各材料溶解率wd,j，所述的各燃料组件各材料溶解率计算如公式(9)：

式中：

dmd,j——b4c、zr、zro2、ss和sso溶解质量；

dt——时间步长；

dmd,j其具体的表达式如公式(10)：

dmd,j＝nmolwj(10)

式中：

n——材料的摩尔数；

molwj——材料的分子量；

3)基于2)得到的各燃料组件各材料溶解率，计算堆芯源节点的熔融物质量流率wf，所述的堆芯源节点的熔融物质量流率计算方程如公式(11)：

wf＝rcdufxpn(11)

式中：

rc——熔化堆芯材料质量密度；

d——熔化堆芯材料稳态膜厚度；

uf——熔化堆芯材料平均速度；

xp——熔化堆芯材料当量直径；

n——源节点中的燃料棒数。

如果熔融物流动为膜状，源节点处圆柱燃料棒上稳态膜厚度d按下式(12)计算

式中：

——源节点处熔化堆芯材料质量；

xpl——源节点处堆芯熔化材料当量直径；

4)利用3)中得到的堆芯源节点的熔融物质量流率求解下一时刻的接受节点的冷凝的熔融物质量mfz，所述的接受节点的冷凝的熔融物质量mfz计算如公式(13)：

mfz＝rcxpnldc/2(13)

式中：

rc——熔化堆芯材料质量密度；

dc——熔融物稳态膜厚度；

l——熔化堆芯材料平均速度；

xp——熔融物稳态膜长度；

n——接受节点中的燃料棒数；

逐渐增长的冷凝外壳的计算与在一半无限体上冷凝一个处于其熔点的液体物质的问题相同；对于这种情况，传导理论产生的平方根增长定律用来计算不断增长的冷凝外壳的瞬时厚度dc，如公式(14)：

式中：

ac——熔融物的热扩散系数；

l——增长常数；

t——熔融物熔点；

增长常数用下面的隐式关系(15)给出：

式中：

b——无因次熔化潜热；

s——修正系数；

l——增长常数；

5)基于4)中得到的接受节点的冷凝的熔融物质量，计算接收节点中待重新分布的熔融物质量mac，所述的接收节点中待重新分布的熔融物质量mac计算方程如公式(16)：

mac＝dtw-max(dtrw,mfz)(16)

式中：

w——熔融物的质量流率；

dt——时间步长；

dtr——当前的时间步长；

mfz——一个步长后接受节点冷凝熔融物质量；

值得注意的是这个质量要加入开始时间步长接受节点已存在的熔融物质量中，将两质量和作为熔融物质量来估算同一时间步长内从接收节点流到其下节点的熔融物质量；

6)利用5)计算得到的接收节点流到其下节点的熔融物质量进行下一步堆芯熔化计算，即重复3)、4)、5)、6)的计算过程，直到计算结束；

步骤3：利用步骤1和步骤2计算得到的各燃耗下的堆芯节点的分布温度、堆芯包壳氧化量及氢气产生质量、包壳应力应变的变化率、各燃料组件各材料溶解率、堆芯源节点的熔融物质量流率、接受节点的冷凝的熔融物质量和接收节点中待重新分布的熔融物质量进行堆芯熔化计算，堆芯熔化计算具体包含以下内容：

1)读取堆芯熔化计算得到的几何区域信息以及堆芯节点熔融物状态参数；

2)利用各燃料组件各材料溶解率、堆芯源节点的熔融物质量流率、接受节点的冷凝的熔融物质量及接收节点中待重新分布的熔融物质量，对熔融物碎片床进行初始化计算，首先利用龙格库塔法计算熔融物迁移到下腔室水池过程中的射流碎裂和颗粒化的过程，所述的熔融物迁移到下腔室水池过程中的射流碎裂和颗粒化的计算过程如公式(17)-(21)：

射流碎裂部分的质量流量wcme计算公式(17)：

wcme＝fentwcmtj(17)

式中：

wcmtj——碎片射流总的质量流量；

fent——连续射流的碎裂分量；

ddj——水中夹带的颗粒直径；

ddj,0——水中夹带的颗粒初始直径；

在水中尺寸范围为1-5mm的碎片颗粒的沉降速度很大，通过下式(19)计算：

式中：

——碎片颗粒的沉降速度；

cd——拖曳系数；

g——重力加速度；

ρdp——水中夹带的颗粒密度；

ρw——水的密度；

ddp——水中夹带的颗粒直径；

夹带后的碎片颗粒温度计算如下式(20)：

式中：

e——夹带系数；

cp,dj——水中夹带的颗粒定压热容；

rdp——水中夹带的颗粒密度；

s——颗粒沉降系数；

tsed——沉降时间；

tdj——射流碎裂时的温度；

3)基于2)得到的熔融物迁移到下腔室水池过程中的射流碎裂和颗粒化的状态参数，计算产生的碎片床高度zlp，所述的产生的碎片床高度计算方程如公式(21)：

式中：

rlp——压力容器半径；

zlp——碎片床高度；

ncrd——crd管道数；

rcrd——crd管道外半径；

nins——仪器管道管道数；

rins——仪器管道外半径；

ash——围板横截面积；

4)利用3)中得到的碎片床高度求解下一时刻的碎片床向压力容器壁面的传热率，所述的碎片床向压力容器壁面的传热率计算如公式(22)：

式中：

tdo——熔融碎片的初始温度；

tso——钢壁面的初始温度；

ks——钢壁面热导率；

kc——熔融物热导率；

rs——钢壁面密度；

rc——熔融物密度；

cc——熔融物厚度；

cs——钢壁面厚度；

钢壁面的平均温度为界面温度和钢壁初始温度的平均值，平均硬壳温度则假定为硬壳熔点和界面温度的算术平均值；

5)基于4)中得到的碎片床向压力容器壁面的传热率，使用封闭法则计算碎片床辐射模型角系数，所述的碎片床辐射模型角系数计算方程如公式(23)：

f12＝f21＝1-(x²/2)((1+4/x²)¹²-1)(23)

式中：

x——等于h/r，为圆盘间的无量纲距离；

f12,f21——碎片床辐射模型角系数；

6)利用5)计算得到的碎片床辐射模型角系数进行下一步堆芯熔化计算，即重复3)、4)、5)、6)的计算过程，直到计算结束；

步骤4：利用步骤3计算得到的熔融物迁移到下腔室水池过程中的射流碎裂和颗粒化的过程、产生的碎片床高度、碎片床向压力容器壁面的传热率和碎片床辐射模型角系数进行压力容器下封头内熔融物构型计算，熔融物构型计算具体包含以下内容：

1)读取堆芯碎片床计算得到的几何区域信息以及堆芯碎片床状态参数；

2)利用产生的碎片床高度、碎片床向压力容器壁面的传热率、碎片床辐射模型角系数，对压力容器下封头熔融物构型进行初始化计算，首先利用龙格库塔法计算熔融物与周围物质的传热量，所述的熔融物与周围物质的传热量计算如公式(24)-(26)：

式中：

qu^ⅱ——平均向上热流密度；

qd^ⅱ——平均向下热流密度；

qs^ⅱ——平均侧向热流密度；

kpx——熔融物热导率；

dt——熔融物过热度；

rlp——下封头半径；

zlp——下封头半径；

ra——雷诺数；

在下封头熔融池内衰变热的作用下，氧化物层内的热量同时向上和向下传递；当压力容器处于长期冷却状态时，下封头内熔融池的上部空间会存在一个较为强烈的湍流自然对流区域，而下部空间则是具有明显热分层的稳定区域；

3)基于2)得到的熔融物与周围物质的传热量，计算熔融物在冷却对流过程中释放到堆芯的蒸汽流量wst，所述的熔融物在冷却对流过程中释放到堆芯的蒸汽流量wst计算方程如公式(27)：

式中：

mw——水的质量；

uw——比内能；

cv,w——定容比热；

qcm——碎片床传递给上部水池的热流量；

ww,in——进入水池的水流量；

hin——进入水池的焓值；

pps——主系统压力；

hfg——蒸汽饱和焓值；

hst——水的饱和焓值；

hw——水的比焓值；

tsat——水的饱和温度；

vfg——饱和水和饱和蒸汽间的比体积差；

mg,ps——在主系统压力下的蒸汽质量；

4)利用2)中得到的熔融物与周围物质的传热量求解下一时刻熔融物颗粒与周围环境间的质量和能量交换，所述的熔融物颗粒与周围环境间的质量和能量交换计算如公式(28)：

式中：

k——熔融物颗粒的热导率系数；

kg——熔融物颗粒窄缝气体的热导率；

kp——熔融物颗粒的热导率；

熔融物颗粒单位体积的热容量计算公式为下式(29)：

rc＝argcg+(1-a)rpcp(29)

式中：

rgcg——熔融物颗粒窄缝气体的热导率；

rpcp——熔融物颗粒的热导率；

a——熔融物颗粒窄缝气体占比系数；

瞬时热流密度q^ⅱ的所有数值结果与下式(30)的误差小于10％：

式中：

k——熔融物颗粒的热导率系数；

tw——冷却剂温度；

to——熔融物颗粒温度；

d——熔融物颗粒厚度；

5)利用4)中得到的熔融物颗粒与周围环境间的质量和能量交换求解下一时刻碎片硬壳的热流量ql^ⅱ，所述的碎片硬壳的热流量ql^ⅱ计算如公式(31)：

式中：

ql^ⅱ——代表在x＝l处从碎片主体到硬壳的对流传热量；

ux——内能的变化率；

qo^ⅱ——代表在x＝o处从碎片主体到硬壳的对流传热量；

mx——硬壳质量；

h——对流传质的焓值；

qv——衰变热体积产热率；

a——碎片熔融池与硬壳间的接触面积；

l——硬壳厚度；

6)利用5)计算得到的碎片硬壳的热流量进行下一步熔融物构型分析计算，即重复3)、4)、5)、6)的计算过程，直到计算结束。

与现有技术相比，本发明具有如下突出特点：

1、利用机械学的模型计算燃料棒的弹性‐塑性变形。各向同性的塑性变形则用hill理论和prandtl‐reuss方程。通过迭代计算得出的收敛平均周向应力，要判定包壳是否破裂。如果包壳破裂，那么包壳变形计算在接下来的分析中就不用进行。

2、对于堆芯熔化过程中的低温共晶作用的描述，采用最新的实验研究数据和理论关系式的基础上，针对堆芯材料之间可能形成的低温共晶作用建立相应的理论分析模型。

3、研究节块间传热时，为了避免解时间‐空间耦合的传热方程，假设了一个准稳态。气体质量被假定为通过节块的流量和时间步长的乘积，且时间步长结束时算得的气体温度被取为出口气体温度以简化计算。

附图说明

图1是压水堆严重事故计算流程示意图。

图2是早期堆芯行为及堆芯燃料组件应力-应变计算流程图。

图3是堆芯熔化计算流程图。

图4是堆芯碎片床特性计算流程图。

图5是压力容器下封头内熔融物构型分析计算流程图。

具体实施方式

下面结合附图和具体实施方式对本发明方法进行详细的说明。

如图1所示，本发明压水堆核电厂严重事故分析方法及技术方法，包括如下步骤：

步骤1：进行严重事故早期堆芯行为特性及燃料组件应力-应变特性计算，计算得到各堆芯节点的分布温度、堆芯包壳氧化量、氢气产生质量以及燃料包壳应力应变的变化率，如图2所示，早期堆芯行为特性及燃料组件应力-应变特性计算具体包含以下内容：

1)读取所有事故相态中堆芯的热工水力参数和响应信息；

2)对计算堆芯进行几何建模：根据计算需求对计算区域进行剖分，给出了一个7×13的堆芯节点划分，有三个非燃料层：两个在底部，代表堆芯支撑板和下部固定板以及下部气体腔室；一个在顶部，代表上部气体腔室和上部固定板；其他的层轴向长度相同，代表堆芯活性区；根据计算需求对计算区域进行剖分，在每一个节点内生成初始状态参数，对每个计算子区域的材料中的各状态进行赋初值；3)基于2)所得到的初始状态参数和各材料的状态信息，利用龙格库塔法进行计算得到堆芯各节点的温度分布tn，所述的堆芯各节点温度分布tn的计算如公式(1)：

式中：

qn——流过堆芯节点的热流量；

wst——水质量流量；

hst——在冷却剂温度为tn，分压为时的焓值；

——冷却剂饱和分压；

4)基于3)所得到堆芯各节点温度分布，计算得到堆芯包壳氧化量所述的堆芯包壳氧化量的计算如公式(2)：

式中：

t——堆芯节点温度即包壳，控制棒或者水棒，燃料包壳的温度；

x——氧化物厚度；

ρzr——锆合金密度；

a——氧化系数；

b——温度修正系数；

r——理想气体常数；

5)基于4)所得到堆芯包壳氧化量，计算包壳氧化过程中的氢气产生质量，所述的氢气产生质量如反应式(5)所示，所述的氢气产生质量计算公式如(6)所示：

zr+2h2o→zro2+2h2+δhzr(5)

式中：

δhzr——每摩尔的锆合金反应产生的热量；

hst——蒸汽入口焓值；

——金属温度下的氢气焓值；

——氢气的质量流量；

wst,rct——蒸汽的质量流量；

qrct——反应产生的化学热；

——锆合金消耗的摩尔速率；

6)基于3)和5)所分别得到的堆芯各节点的温度分布和堆芯包壳氧化量，计算得到各堆芯节点的燃料包壳应力应变的变化率，所述的燃料包壳应力应变变化率的计算如公式(7)、(8)：

式中：

σr，σθ，σz——径向应力，周向应力，轴向应力；

εr，εθ，εz——径向应变，周向应变，轴向应变；

——径向塑性应变量，周向塑性应变量，轴向塑性应变量；

——径向塑性应变增量，周向塑性应变增量，轴向塑性应变增量；

——径向热变形应变量，周向热变形应变量，轴向热变形应变量；

——径向弹性变形应变量，周向弹性变形应变量，轴向弹性变形应变量；

e——杨氏模量；

v——泊松比；

αr，αθ，αz——径向热膨胀率，周向热膨胀率，轴向热膨胀率；

——膨胀导致的热应变；

boz——boltzmann常数；

k1,k2,k3,ed——经验常数；

t——温度；

由于燃料棒的长度远大于其半径尺寸，计算时将其视为一维平面应力问题进行处理；

7)利用6)计算得到各堆芯节点的燃料包壳应力应变的变化率进行下一步早期堆芯行为特性及燃料组件应力-应变特性计算，即重复3)、4)、5)、6)的计算过程，直到计算结束；

步骤2：基于步骤1中计算得到的各堆芯节点下的温度分布、堆芯包壳氧化量、氢气产生质量以及燃料包壳应力应变的变化率进行堆芯熔化计算，堆芯熔化计算相比于早期堆芯行为特性及燃料组件应力-应变特性计算在时间上是一个顺向的过程，堆芯熔化计算具体包含以下内容：

1)读取早期堆芯行为特性及燃料组件应力-应变特性计算得到的几何区域信息以及堆芯节点各状态参数；

2)利用各堆芯节点下的温度分布、堆芯包壳氧化量、氢气产生质量以及燃料包壳应力应变的变化率，对燃料组件各参数进行初始化计算，首先计算各燃料组件各材料溶解率wd,j，如图3所示，所述的各燃料组件各材料溶解率计算如公式(9)：

式中：

dmd,j——b4c、zr、zro2、ss和sso溶解质量；

dt——时间步长；

dmd,j其具体的表达式如公式(10)：

dmd,j＝nmolwj(10)

式中：

n——材料的摩尔数；

molwj——材料的分子量；

3)基于2)得到的各燃料组件各材料溶解率，计算堆芯源节点的熔融物质量流率wf，所述的堆芯源节点的熔融物质量流率计算方程如公式(11)：

wf＝rcdufxpn(11)

式中：

rc——熔化堆芯材料质量密度；

d——熔化堆芯材料稳态膜厚度；

uf——熔化堆芯材料平均速度；

xp——熔化堆芯材料当量直径；

n——源节点中的燃料棒数。

如果熔融物流动为膜状，源节点处圆柱燃料棒上稳态膜厚度d按下式(12)计算

式中：

——源节点处熔化堆芯材料质量；

xpl——源节点处堆芯熔化材料当量直径；

4)利用3)中得到的堆芯源节点的熔融物质量流率求解下一时刻的接受节点的冷凝的熔融物质量mfz，所述的接受节点的冷凝的熔融物质量mfz计算如公式(13)：

mfz＝rcxpnldc/2(13)

式中：

rc——熔化堆芯材料质量密度；

dc——熔融物稳态膜厚度；

l——熔化堆芯材料平均速度；

xp——熔融物稳态膜长度；

n——接受节点中的燃料棒数；

逐渐增长的冷凝外壳的计算与在一半无限体上冷凝一个处于其熔点的液体物质的问题相同；对于这种情况，传导理论产生的平方根增长定律用来计算不断增长的冷凝外壳的瞬时厚度dc，如公式(14)：

式中：

ac——熔融物的热扩散系数；

l——增长常数；

t——熔融物熔点；

增长常数用下面的隐式关系(15)给出：

式中：

b——无因次熔化潜热；

s——修正系数；

l——增长常数；

5)基于4)中得到的接受节点的冷凝的熔融物质量，计算接收节点中待重新分布的熔融物质量mac，所述的接收节点中待重新分布的熔融物质量mac计算方程如公式(16)：

mac＝dtw-max(dtrw,mfz)(16)

式中：

w——熔融物的质量流率；

dt——时间步长；

dtr——当前的时间步长；

mfz——一个步长后接受节点冷凝熔融物质量；

值得注意的是这个质量要加入开始时间步长接受节点已存在的熔融物质量中，将两质量和作为熔融物质量来估算同一时间步长内从接收节点流到其下节点的熔融物质量；

6)利用5)计算得到的接收节点流到其下节点的熔融物质量进行下一步堆芯熔化计算，即重复3)、4)、5)、6)的计算过程，直到计算结束；

步骤3：利用步骤1和步骤2计算得到的各燃耗下的堆芯节点的分布温度、堆芯包壳氧化量及氢气产生质量、包壳应力应变的变化率、各燃料组件各材料溶解率、堆芯源节点的熔融物质量流率、接受节点的冷凝的熔融物质量和接收节点中待重新分布的熔融物质量进行堆芯熔化计算，堆芯熔化计算具体包含以下内容：

1)读取堆芯熔化计算得到的几何区域信息以及堆芯节点熔融物状态参数；

2)利用各燃料组件各材料溶解率、堆芯源节点的熔融物质量流率、接受节点的冷凝的熔融物质量及接收节点中待重新分布的熔融物质量，对熔融物碎片床进行初始化计算，首先利用龙格库塔法计算熔融物迁移到下腔室水池过程中的射流碎裂和颗粒化的过程，如图4所示，所述的熔融物迁移到下腔室水池过程中的射流碎裂和颗粒化的计算过程如公式(17)-(21)：

射流碎裂部分的质量流量wcme计算公式(17)：

wcme＝fentwcmtj(17)

式中：

wcmtj——碎片射流总的质量流量；

fent——连续射流的碎裂分量；

ddj——水中夹带的颗粒直径；

ddj,0——水中夹带的颗粒初始直径；

在水中尺寸范围为1-5mm的碎片颗粒的沉降速度很大，通过下式(19)计算：

式中：

——碎片颗粒的沉降速度；

cd——拖曳系数；

g——重力加速度；

ρdp——水中夹带的颗粒密度；

ρw——水的密度；

ddp——水中夹带的颗粒直径；

夹带后的碎片颗粒温度计算如下式(20)：

式中：

e——夹带系数；

cp,dj——水中夹带的颗粒定压热容；

rdp——水中夹带的颗粒密度；

s——颗粒沉降系数；

tsed——沉降时间；

tdj——射流碎裂时的温度；

3)基于2)得到的熔融物迁移到下腔室水池过程中的射流碎裂和颗粒化的状态参数，计算产生的碎片床高度zlp，所述的产生的碎片床高度计算方程如公式(21)：

式中：

rlp——压力容器半径；

zlp——碎片床高度；

ncrd——crd管道数；

rcrd——crd管道外半径；

nins——仪器管道管道数；

rins——仪器管道外半径；

ash——围板横截面积；

4)利用3)中得到的碎片床高度求解下一时刻的碎片床向压力容器壁面的传热率，所述的碎片床向压力容器壁面的传热率计算如公式(22)：

式中：

tdo——熔融碎片的初始温度；

tso——钢壁面的初始温度；

ks——钢壁面热导率；

kc——熔融物热导率；

rs——钢壁面密度；

rc——熔融物密度；

cc——熔融物厚度；

cs——钢壁面厚度；

钢壁面的平均温度为界面温度和钢壁初始温度的平均值，平均硬壳温度则假定为硬壳熔点和界面温度的算术平均值；

5)基于4)中得到的碎片床向压力容器壁面的传热率，使用封闭法则计算碎片床辐射模型角系数，所述的碎片床辐射模型角系数计算方程如公式(23)：

f12＝f21＝1-(x²/2)((1+4/x²)^1/2-1)(23)

式中：

x——等于h/r，为圆盘间的无量纲距离；

f12,f21——碎片床辐射模型角系数；

6)利用5)计算得到的碎片床辐射模型角系数进行下一步堆芯熔化计算，即重复3)、4)、5)、6)的计算过程，直到计算结束；

步骤4：利用步骤3计算得到的熔融物迁移到下腔室水池过程中的射流碎裂和颗粒化的过程、产生的碎片床高度、碎片床向压力容器壁面的传热率和碎片床辐射模型角系数进行压力容器下封头内熔融物构型计算，如图5所示，熔融物构型计算具体包含以下内容：

1)读取堆芯碎片床计算得到的几何区域信息以及堆芯碎片床状态参数；

2)利用产生的碎片床高度、碎片床向压力容器壁面的传热率、碎片床辐射模型角系数，对压力容器下封头熔融物构型进行初始化计算，首先利用龙格库塔法计算熔融物与周围物质的传热量，所述的熔融物与周围物质的传热量计算如公式(24)-(26)：

式中：

qu^ⅱ——平均向上热流密度；

qd^ⅱ——平均向下热流密度；

qs^ⅱ——平均侧向热流密度；

kpx——熔融物热导率；

dt——熔融物过热度；

rlp——下封头半径；

zlp——下封头半径；

ra——雷诺数；

在下封头熔融池内衰变热的作用下，氧化物层内的热量同时向上和向下传递；当压力容器处于长期冷却状态时，下封头内熔融池的上部空间会存在一个较为强烈的湍流自然对流区域，而下部空间则是具有明显热分层的稳定区域；

3)基于2)得到的熔融物与周围物质的传热量，计算熔融物在冷却对流过程中释放到堆芯的蒸汽流量wst，所述的熔融物在冷却对流过程中释放到堆芯的蒸汽流量wst计算方程如公式(27)：

式中：

mw——水的质量；

uw——比内能；

cv,w——定容比热；

qcm——碎片床传递给上部水池的热流量；

ww,in——进入水池的水流量；

hin——进入水池的焓值；

pps——主系统压力；

hfg——蒸汽饱和焓值；

hst——水的饱和焓值；

hw——水的比焓值；

tsat——水的饱和温度；

vfg——饱和水和饱和蒸汽间的比体积差；

mg,ps——在主系统压力下的蒸汽质量；

4)利用2)中得到的熔融物与周围物质的传热量求解下一时刻熔融物颗粒与周围环境间的质量和能量交换，所述的熔融物颗粒与周围环境间的质量和能量交换计算如公式(28)：

式中：

k——熔融物颗粒的热导率系数；

kg——熔融物颗粒窄缝气体的热导率；

kp——熔融物颗粒的热导率；

熔融物颗粒单位体积的热容量计算公式为下式(29)：

rc＝argcg+(1-a)rpcp(29)

式中：

rgcg——熔融物颗粒窄缝气体的热导率；

rpcp——熔融物颗粒的热导率；

a——熔融物颗粒窄缝气体占比系数；

瞬时热流密度q^ⅱ的所有数值结果与下式(30)的误差小于10％：

式中：

k——熔融物颗粒的热导率系数；

tw——冷却剂温度；

to——熔融物颗粒温度；

d——熔融物颗粒厚度；

5)利用4)中得到的熔融物颗粒与周围环境间的质量和能量交换求解下一时刻碎片硬壳的热流量ql^ⅱ，所述的碎片硬壳的热流量ql^ⅱ计算如公式(31)：

式中：

ql^ⅱ——代表在x＝l处从碎片主体到硬壳的对流传热量；

ux——内能的变化率；

qo^ⅱ——代表在x＝o处从碎片主体到硬壳的对流传热量；

mx——硬壳质量；

h——对流传质的焓值；

qv——衰变热体积产热率；

a——碎片熔融池与硬壳间的接触面积；

l——硬壳厚度；

6)利用5)计算得到的碎片硬壳的热流量进行下一步熔融物构型分析计算，即重复3)、4)、5)、6)的计算过程，直到计算结束。

大量计算验证结果显示，本发明具有可靠的精度、很好的效率和很好的几何适应性，适应工程实际中的计算要求。程序可以从电厂正常运行开始，计算瞬态运行，loca事故，严重事故序列分析(早期行为，堆芯融化，碎片床，堆内熔融物保持分析)。程序编制考虑了更多的物理模型，并采用了更贴近实际堆型的经验关系式，比较了三大严重事故分析程序的模型特点，添加了更为准确的数值计算方法，提高了计算的精度，所得的计算结果更加可靠。