1. 引言
《核动力厂设计安全规定》(HAF102)对严重事故的预防与缓解有明确要求。一旦发生严重事故，将会产生大量的堆芯熔融物，这些熔融物碎片不仅具有相当强的放射性，而且还会继续衰变产生衰变热。因此，研究堆芯熔融物在一回路内的迁移及分布特性, 对于准确计算安全壳源项分布，研究严重事故缓解措施及应急策略具有重要意义。
熔融物碎片在冷却剂内的流动特性是固液两相流问题。由于熔融物碎片尺寸较大，传统的网格法在求解此类问题时容易导致网格变形，且商用计算软件中的VOF模型、欧拉模型、DPM模型等也不再适用。无网格方法由于其追踪粒子方便，不依赖于网格而被广泛地应用到天体物理学及形变剧烈的问题中 [1] 。1996年，Koshizuka和Oka等人提出了移动粒子半隐式方法(MPS) [2] [3] ，该方法可以准确地跟踪离散流体的运动及捕捉自由液面，并扩展到了多相多流体的各种热工水力学的问题中，比如自由表面的流动 [4] 、液滴的破裂 [5] 、核态沸腾时气泡的产生、脱离以及上升 [6] 与在高热流密度和过冷度时的沸腾 [7] 、射流的喷射 [8] 、气体的破裂 [9] 和液体的破裂以及固体发生大变形与破裂 [10] 等现象。西安交通大学的田文喜、陈荣华、苏光辉等 [11] [12] [13] 应用MPS-MAFL方法研究氩气在液态铅铋合金内上升过程中的二维运动特性，并考虑了流体间的换热、气泡的形状等因素,同时基于传统的两相流理论研究了气泡提升泵工作回路内空泡份额对自然循环能力的影响 [14] 。华北电力也开展了基于粒子法的快堆主容器晃动问题的研究 [15] 。
本文在粒子法程序MPS中增加铅/铅铋流体计算模型，并定义其产生条件、入口速度、重力、粘性力、粒子数密度计算模型。利用改进后的程序模拟Zr-4合金碎片及二氧化铀碎片在不同的冷却剂流速及冷却剂粘性条件下的运动行为，初步验证粒子法程序的有效性，为研究铅冷快堆在严重事故工况下熔融碎片的流动及分布情况提供了程序基础。
2. MPS理论与方法
MPS方法基于拉格朗日方法对流场进行描述，通过定义粒子间的相互作用用来描述梯度，散度，拉普拉斯算子等模型，从而对N-S方程进行离散求解。在MPS方法中，采用以下核函数来表示粒子间的相互作用(图1)。
(1)
式中：r为相邻两粒子间距，r越小，相互作用越大；re为核尺寸，决定了粒子间相互作用的面积大小。
粒子i和粒子j的位置用ri，rj来表示。以粒子i为中心，以i和j之间的矢量差r = |rj − ri|为变量，对它进行加权求和得到了粒子数密度。
(2)
对于不可压流体，流体密度不变相当于粒子数密度不变。将梯度、散度、拉普拉斯算子用粒子间的相互作用表示，就可以对连续性方程以及N-S方程进行离散。
粒子i处的梯度矢量计算模型为：
(3)
粒子i (d维空间)的散度公式如下：
(4)
通过拉普拉斯算子表示控制方程中物理量的空间扩散项，其表达式位：
(5)
MPS粒子法采用半隐式的时间推进算法来计算不可压缩流体的流动。首先显式计算出颗粒的位置和速度，然后求解控制方程，得到速度和位置的修正值，最后经过修正第一步的结果，得到下一时刻各个参数的准确值。
对流项的稳定性由Courant数来控制，扩散项的稳定性由扩散数来控制。Courant数的定义如下：
(6)
其中，ui 代表粒子i速度的绝对值，l0 代表粒子初始间距。为了保证数值计算的稳定性，要求：
(7)
3. MPS方法对于固液两相流模拟的有效性验证
在MPS粒子法的半隐式算法中，压力梯度项为
，其与流体密度成反比，因此，如果两相流中的两种流体其密度比过大，界面上很小的压力变化就会给密度小的流体粒子施加巨大的加速度变化，
Figure 1. Interaction between particles
图1. 粒子间相互作用图
从而使其剧烈运动，导致界面上计算的不稳定。为了避免计算不稳定，两种流体的密度比不宜过大。根据以往计算经验，当两种流体密度比在10以内时，程序可以进行稳定的计算，本文中熔融碎片与水冷却剂的密度比值为9，熔融碎片与铅冷却剂的密度比值约为1.16。
严重事故工况核反应堆堆芯熔融碎片的主要成分有核燃料二氧化铀的碎片、燃料包壳熔融再凝固形成的Zr合金碎片等，在一回路温度下，
的密度在10,000 kg/m3左右，而Zr-4合金的密度在6000 kg/m3左右，本计算取堆芯熔融碎片的密度为9000 kg/m3。
在MPS程序中加入了流体种类代号为4的熔融物颗粒，并定义其产生条件、入口速度、重力(cal_buoyancy1函数)、对流项(cal_convection1函数)、粒子数密度(cal_fldnumb函数)计算模型。当判断流体种类为4，即为图2中红色流体粒子时，为其赋值一个x方向的速度，这样程序就能够定义入口处的流体流速。对流项函数用于计算4号粒子的对流运动，如公式(8) (9)所示，其中，
就是上文中定义的入口处流速。
(8)
(9)
粒子数密度的计算公式相同，但是包含的粒子数种类不同，如式(10)所示：
(10)
应用MPS方法建立的计算模型如图2所示，管道的壁面由绿色与黑色粒子表示，其中绿色粒子代表管道内壁，在计算时具有压力，外侧的黑色粒子代表管道外壁，计算时不考虑其压力，由于本模拟中权函数的半径 取值为4.0倍的粒子间距离，故布置两层黑色粒子。粉色粒子代表冷却剂流体，在模拟工况1中为密度为1000 kg/m3的水，而在模拟工况2中，其代表密度为10,419 kg/m3的液态铅冷却剂。红色粒子代表熔融物碎片与冷却剂的混合物，布置在管道入口处。熔融物碎片进入管道内部后会变成蓝色粒子密度为9000 kg/m3。流体的入口流速参考压水堆核电站冷却剂平均流速，设置为4.65 m/s。
管道模型的出口边界为压力为0的开放性出口，入口边界为速度入口，程序运行时新的流体粒子不断从入口处流入，推动管道内的粒子从右侧出口流出，这样便模拟了颗粒物在水中的流动现象。
通过对比图3工况(a)与工况(b)的结果可得看出，当两种粒子密度差较大时，固体粒子比液体粒子更容易沉积在管道底部和障碍物处。而当密度差不大时，固体粒子会与液体粒子混杂在一起沿着管道运动，并且不会发生明显的沉积现象。该结果与实际经验相符，初步证明了MPS程序在模拟固液两相流现象上的可行性。
4. 不同种类的熔融碎片在铅冷却剂内的流动模拟
建立与图3相同的计算模型，设置管道内径为0.7米，壁面粒子为316钢，密度为7850 kg/m3，参考550℃下铅的密度为10,419 kg/m3。选取二氧化铀、Zr-4合金两种熔融物碎片，模拟不同熔融物碎片、粘度(温度)、流速等条件下碎片的运动情况。
4.1. 不同熔融物碎片的流动行为
熔融碎片为Zr-4合金时，将固体粒子密度设为6000 kg/m3；熔融碎片为二氧化铀时，将固体粒子密度设为10,000 kg/m3。
图4显示的是Zr合金碎片与二氧化铀碎片在铅冷却剂内的迁移行为模拟图。从图4可以看出，当碎片成分是Zr合金时，由于其密度较小，在重力的作用下，流动时会逐渐漂浮到冷却剂流体的上面。当碎片成分为二氧化铀时，由于其密度与冷却剂流体的密度十分相近，加上本文采用的固体粒子与液体粒子同等尺寸，因此二氧化铀粒子与流体粒子混杂在一起，难以发生分离和沉积现象，并且即使当流速改变的时候其运动状况也没有明显的变化。
(a) 混合熔融物碎片在水内的流动
(b) 混合熔融物碎片在铅内的流动
Figure 3. Simulation of molten debris in coolant with different densities
图3. 混合熔融物颗粒在不同密度冷却剂内的流动模拟
4.2. 不同冷却剂流速条件下熔融物碎片的流动行为
图5显示的是冷却剂流速分别为0.5 m/s、1 m/s、2 m/s时Zr合金碎片在铅冷却剂内的运动行为。当冷却剂流速较小时，Zr合金碎片在重力的作用下会逐渐漂浮到冷却剂流体的上层。随着冷却剂流速的增加，Zr合金碎片上浮的现象逐渐减弱，而更倾向于保持原有的运动状态。这是因为流速越大，冷却剂粒子对固体碎片的夹带作用越明显。
4.3 不同冷却剂运动粘度条件下熔融物碎片的流动行为
事故工况下冷却剂的温度将发生较大的波动，会导致冷却剂的粘性发生变化，因此考虑冷却剂有粘性和无粘性两种工况，计算二氧化铀碎片和Zr合金碎片在冷却剂内的运动行为，如图6和图7所示。
当冷却剂流速为0.5 m/s且不考虑冷却剂粘性时，虽然二氧化铀碎片(图6-a)和Zr碎片(图7-a)存在密
(a) 流速为0.5 m/s时Zr合金的流动
(b) 流速为0.5 m/s时二氧化铀碎片的流动
Figure 4. Simulation of molten debris in coolant with different densities
图4. 混合熔融物颗粒在不同密度冷却剂内的流动模拟
(a) 流速为0.5 m/s时二氧化铀碎片的流动(无粘)
(b) 流速为0.5 m/s时二氧化铀碎片的流动(有粘)
Figure 6. Simulation of UO2 debris in lead coolant with different viscosities
图6. 二氧化铀碎片在铅冷却剂内的流动
(a) 流速为0.5 m/s时Zr合金的流动(无粘)
(b) 流速为0.5 m/s时Zr合金的流动(有粘)
Figure 7. Simulation of Zr-4 alloy debris in lead coolant with different viscosities
图7. Zr合金碎片在铅冷却剂内的流动
度上的差别，但两种碎片与冷却剂都有不同程度的搅混。考虑冷却剂粘性以后，两种碎片与冷却剂的搅混程度明显减弱，如图6-b和图7-b所示。这是由于冷却剂粘性导致其对固体粒子的夹带作用明显增加，并且使固体粒子发生搅混的阻力增大。
5. 结论
本文通过定义熔融物颗粒的产生条件、入口速度、重力、粘性力、粒子数密度计算模型，完成粒子法程序MPS中铅/铅铋流体计算模型的开发。利用改进后的程序计算Zr-4合金碎片及二氧化铀碎片在不同的冷却剂流速及冷却剂粘性条件下的运动行为。
冷却剂流速较低时，Zr-4合金碎片由于其密度较小，在流动时会逐渐上浮到冷却剂的上层。而二氧化铀颗粒与流体粒子混杂在一起，难以发生分离和沉积现象。随着冷却剂流速的增加，冷却剂粒子对固体碎片的夹带作用明显，合金碎片上浮的现象逐渐减弱，而更倾向于保持原有的运动状态。考虑冷却剂粘性后，固体粒子发生搅混的阻力增大，两种碎片随冷却剂运动的一致性程度增加。
作为应用粒子法开展严重事故工况下熔融物运动行为研究的初步适用性验证阶段，本文仍采用固体粒子与液体粒子等粒径模型，下一阶段将开发不同尺寸的熔融物碎片运动行为计算模型，并开展相关的实验验证。
基金项目
国际原子能机构资助项目(R18971)。