曹 胜,张 斌,2,*,王文鹏,单建强,2
(1.西安交通大学 能源与动力工程学院,陕西 西安 710049;2.西安交通大学 动力工程多相流国家重点实验室,陕西 西安 710049)
钠冷快堆(SFR)堆芯解体事故(CDA)可能对堆内构件和压力容器造成严重损害,是放射性包容与纵深防御的重要安全隐患。当熔融物下降并与下腔室中较低温度的冷却剂接触后,会骤冷形成放射性碎片,最终沉降在反应堆下腔室形成碎片床[1]。若碎片床没有被妥善地滞留和冷却,可能导致压力容器失效。因此,钠冷快堆中设计了堆芯捕集器作为预防与缓解装置,用于收集和容纳熔融物碎片,以尽可能地实现对碎片床的长期冷却,确保主容器的完整性。
目前已有不少学者研究了碎片床的冷却特性和堆芯捕集器的性能。Sudha等[2]通过数值传热分析估算了钠冷快堆在无保护失流事故(ULOFA)和热阱丧失事故(PLOHS)下熔融物转移到堆芯捕集器的最短和最长时间,这些估算时间对于堆芯捕集器的性能设计具有重要的参考意义。David等[3]分析了SFR池内自然对流的瞬态发展,并评估了不同堆芯熔化场景下堆芯捕集器底板的最高温度。Roychowdhury等[4]使用一维模型证明了印度原型快堆(PFBR)能够安全消除19个组件融化带来的放射性衰变热。Vlasichev等[5]使用STO-BEO代码对BN800捕集器上的碎片床开展了数值分析,并评估出碎片床内的最高温度低于燃料沸腾温度。在许多研究中碎片床的形状是任意选择的,如圆柱形、圆锥形、高斯形、堆状、圆锥状或同向堆积[6]。然而,碎片床的形状和分布很重要,它是防止碎片床重返临界和衰变热长期冷却的关键影响因素。
一些快堆在设计和建造阶段就依据反应堆特征对堆芯捕集器进行了优化。例如,英国Dounreay快堆的下腔室中设计了燃料分散锥和防溅板,以避免堆芯熔融物结块,并将熔融物引导至堆芯捕集器[7-8];英国商业示范快堆(CDFR)采用多个垂直管道焊接成圆锥形管板收集堆芯碎片,并连接到一个中央烟囱以促进冷却剂的自然对流[8-9];俄罗斯BN800和BN1200所配备的堆芯捕集器,采用7个垂直引流管以增强自然对流[10];欧洲快堆(EPR)配备了箱式结构的堆芯捕集器,并采用双烟囱结构以更好地分散堆芯碎片[9,11]。此外,还有采用多托盘的堆芯捕集器,将堆芯碎片分开进行冷却[12-14],采用多层托盘的堆芯捕集器设计,消除了冷却剂对堆芯捕集器部件长期侵蚀和腐蚀的隐患[7,15]。这些设计通过工程实践对堆芯捕集器的结构及组件进行了优化。烟囱作为堆芯捕集器的主要结构之一,不仅影响冷却剂的自然对流,而且对碎片床的堆积和分布具有重要影响,因此亟需进行详细的研究与设计优化。
针对目前堆芯捕集器烟囱结构设计优化研究不足的现状,本文采用改进离散元法(DEM)作为数值模拟工具,以中国实验快堆(CEFR)为原型,对碎片床在堆芯捕集器上的形成过程进行数值模拟。通过改变堆芯捕集器托盘上的烟囱顶盖垂直投影边长、烟囱顶盖倾角以及烟囱间距,研究堆芯捕集器烟囱结构对碎片床堆积形状和分布的作用机理和影响规律,为碎片床的传热性能分析和长期可冷却性评估提供前提条件,并对钠冷快堆堆芯捕集器的设计及优化提供工程参考。
本文采用DEM模拟堆芯捕集器上碎片床的形成过程。在DEM计算中,将熔融物碎片假设为粒子,其运动遵循牛顿定律。碎片颗粒的运动分为平移和旋转,并通过以下控制方程分别进行求解:
(1)
(2)
式中:mi和Ii为颗粒i的质量和转动惯量;ui和ωi为颗粒i的平移速度和旋转速度。两粒子碰撞时粒子i的受力如图1所示,其中Fi和Ti分别为作用在粒子i上的力和力矩,由方程(3)和(4)求解。
(3)
图1 粒子i与粒子j碰撞的受力示意图
(4)
粒子间和粒子与壁面间碰撞产生的法向接触力Fn和切向接触力Fs采用弹簧-阻尼模型计算,如式(5)、(6)所示,图2为弹簧-阻尼模型示意图。
(5)
(6)
式中:δ为粒子接触点的相对位移;k为刚度系数;c为阻尼系数;μ为摩擦系数。
法向刚度系数kn和阻尼系数cn分别由式(7)和式(8)计算。
(7)
(8)
式中:αtn为无量纲刚度修正系数;αcn为无量纲阻尼修正系数;Δt为数值模拟中的时间步长。
数值模拟中,首先确定无量纲刚度系数、无量纲阻尼系数和时间步长的数值。在DEM中,时间步长的选取需要考虑计算效率和计算精度,并结合计算机的运行性能满足公式Δt≤2(m/kn)1/2。通过单个颗粒自由落体后的回弹高度确定无量纲刚度系数αtn的值。通过计算1组颗粒沉降算例,限制粒子堆积的重叠度小于颗粒直径的1%,确定无量纲阻尼系数αcn[16]。使用αtn、αcn计算kn和cn后,通过式(9)、(10)计算切向刚度系数ks和阻尼系数cs。
(9)
(10)
式中,υ为泊松比。在均匀等向性材料中,剪切模量G、杨氏模量E与υ满足G=E/2(1+υ)。
本文采用αtn和αcn替代传统DEM计算中的kn、ks、cn和cs,不仅减少了需要确定的参数数量,也减少了确定参数取值带来的工作量。
除颗粒间的碰撞外,颗粒也受到冷却剂的曳力Fdrag、浮力Fbuoy以及颗粒自身的重力Fg,分别通过式(11)、(12)和(13)计算[17]。
(11)
(12)
(13)
式中:d为颗粒的等效直径;ρd为颗粒的密度;ρc为冷却剂的密度;u为颗粒与冷却剂的相对速度;g为重力加速度;Cdrag为曳力系数,是绕流雷诺数Re的函数。
本文采用编程方式进行DEM求解,计算流程如图3所示。在DEM计算中,采用显式方案更新粒子的速度和位置,计算数据通过开源软件Tecplot进行后处理,以可视化颗粒的运动和碎片床的形成过程。
本文以CEFR堆芯捕集器为研究对象。在CEFR堆芯捕集器设计中,考虑了两种事故假设:1) 假设转化区和活性区的所有燃料组件熔化并迁移至堆芯捕集器;2) 假设1/2的燃料组件迁移到堆芯捕集器。事故参数如下:活性区燃料质量为13.439 t,转化区燃料质量为13.093 t,包壳与盒壁的质量为13.894 t。假设事故后混合的熔融物体积不变,约为4.3 m3。为确保安全裕度,假设事故形成的碎片床孔隙率为0.5,堆芯捕集器的托盘设计容积为9.00 m3。综合考虑数值模拟的计算效率和验证实验的对照等因素后,本研究采用比例为1∶8的模型对堆芯捕集器进行建模和计算。图4为堆芯捕集器的建模图纸,该装置由带有5个烟囱的底部托盘和正上方的“漏斗”状释放装置组成,烟囱设置为1个中央烟囱和4个周边烟囱,且周边烟囱呈等距排布。考虑到堆芯捕集器具有旋转对称性,仅对堆芯捕集器的1/4进行建模和计算,建模示意如图5所示。堆芯捕集器结构由黑色粒子组成,熔融物碎片颗粒由绿色粒子构成,粒子直径为6.0 mm[18]。熔融物粒子静止于释放装置内,在模拟开始后依靠重力自由下落。数值模拟中熔融物粒子和堆芯捕集器壁面粒子的参数如下:密度均为3 700.00 kg/m3,泊松比均为0.27,αtn均为300,αcn均为0.041,静摩擦系数均为0.30。
a——整体结构;b——底部托盘俯视图
图5 堆芯捕集器建模
在钠冷快堆严重事故中,冷却剂通过阻力和浮力影响碎片的下落行为,本研究仅考虑冷却剂对堆芯碎片的单向作用。通过在碎片颗粒的动量方程中加入冷却剂的曳力和浮力,考虑冷却剂对堆芯碎片的影响。曳力计算详见式(11),其中曳力系数计算公式如式(14)所示。
(14)
式中:ud和uc分别为碎片颗粒和钠冷却剂的速度;cdrag为阻力系数,与绕流雷诺数Re有关,通过式(15)计算。在本研究中,钠冷却剂视为大空间内流体,冷却剂密度为ρc,碎片颗粒的等效直径为d。
(15)
(16)
由于堆芯碎片的密度远大于冷却剂的密度,因此认为碎片的运动处于平方阻力区,取k=0.48、n=0[19]。钠冷却剂的物性参数如下:液位为0.34 m、密度为771.8 kg/m3、动力黏度为0.001 005 N·s/m2。
本研究通过开展数值模拟验证实验,以证明所使用数值方法具备堆芯捕集器设计与优化的研究能力。图6为验证实验装置示意图,该装置由装置主体、辅助支撑台架、数据采集系统和进排水系统组成。装置主体由透明有机玻璃制成的底部托盘和不锈钢漏斗释放装置组成,如图7所示。底部烟囱由透明有机玻璃制成,为方便更换,与底部托盘采用分离式设计,通过磁铁和固定销固定。数据采集系统由两套装置组成,一套是高速摄像机,用于连续拍摄碎片的下落和迁移过程;另一套是三维激光扫描仪,在颗粒堆积稳定后对碎片床进行扫描,获取相关的堆积形状数据。
图6 堆芯捕集器实验平台示意图
图7 堆芯捕集器的装置主体
实验使用直径6 mm的氧化铝球形颗粒代替碎片颗粒,总体积为7 L。堆芯捕集器底部托盘内的中央烟囱与周围烟囱的距离设定为120 mm,烟囱顶盖的垂直投影边长为60 mm,顶盖倾角为15°,如图8所示。实验开始前,将颗粒倒入释放装置中并抚平表面,向底部托盘内注入指定高度的水。待水面静止后,迅速释放漏斗内的颗粒。待颗粒堆积稳定后,通过排水口缓慢排出托盘内的水。使用三维激光扫描仪对托盘内形成的碎片床进行扫描,得到碎片床堆积形状和颗粒分布数据。实验过程中高速摄像机拍摄到的照片如图9所示,三维激光扫描仪扫描得到的颗粒分布如图10所示,其中缺失部分是烟囱占据的空间。
图8 验证实验中烟囱的结构与参数
图9 高速摄像机拍摄的不同时刻的粒子运动
图10 碎片床稳定堆积后的分布
为定量比较数值模拟结果与实验结果之间的差异,选取3个指标评价碎片床颗粒分布的均匀性,分别是碎片床的高度方差Zvar、高度归一化方差Znorvar和三维模化方差σ2,计算公式如式(17)~(19)所示,它们可以反映碎片床在高度方向和空间维度上的不均匀程度。
(17)
(18)
(19)
将数值模拟获得的碎片床表面粒子坐标和三维激光扫描仪扫描的粒子坐标按照式(17)、(19)计算,结果列于表1。从表1可看到,数值模拟结果与实验结果之间的相对误差最大值为15.0%,在可接受范围内,表明本研究所采用的数值模拟方法能够有效模拟堆芯捕集器上碎片床的形成过程。
表1 数值模拟和实验中颗粒床粒子群的高度方差和三维模化方差
在数值模拟研究中,通过改变烟囱顶盖垂直投影边长、顶盖倾角和烟囱间距探究堆芯捕集器的烟囱设计对碎片床形状和分布的作用机理及影响规律,从而为堆芯捕集器的结构设计和性能优化提供建议和指导。
在数值模拟中,保持中央烟囱与周围烟囱间距为120 mm、烟囱顶盖倾角为15°,通过改变烟囱顶盖的垂直投影边长(50~100 mm)形成第1组模拟工况(工况1a~1g),如图11所示。在模拟完成后,得到碎片床粒子的高度方差和高度归一化方差,结果如图12所示。从图12可观察到,随着烟囱顶盖垂直投影边长的增加,碎片床的高度方差和高度归一化方差呈相同的变化趋势。当烟囱顶盖垂直投影边长为70 mm时,碎片床的高度方差和高度归一化方差达到最小,分别为1 609 mm2和5.394×10-4。这表明在该组工况中,烟囱顶盖垂直投影边长为70 mm时,堆芯捕集器表现出最佳的碎片颗粒分散性能。
图11 烟囱结构参数示意图
图12 碎片床的高度方差和高度归一化方差
对比发现,碎片颗粒在不同的烟囱顶盖上运动及后续行为具有差异。在工况1a中,没有观察到颗粒的二次散射,颗粒从烟囱顶盖滑落后于烟囱底部附近聚集。而在工况1g中,由于烟囱顶盖过大导致烟囱间缝隙较小,使得颗粒下落时在顶盖边缘处滞留并在烟囱缝隙下方发生明显堆积。在工况1c中,烟囱底部附近没有出现碎片床的凹陷,烟囱缝隙下的颗粒堆积也不明显,烟囱顶盖使部分颗粒在顶盖上发生二次散射,弹射到更远的位置,从而形成最均匀的碎片床分布。图13为工况1a~1g中碎片床在y-z平面上的投影。由图13可观察到,随着烟囱顶盖垂直投影边长的增加,滞留在堆芯捕集器烟囱顶盖上的碎片颗粒数量也增加。碎片颗粒滞留的原因是最后下落的颗粒不再受后续颗粒的撞击,导致摩擦力大于沿烟囱切线方向的重力分量。所有工况中碎片床的堆积角均约为32°,但碎片颗粒分布直径D(D=2((x2+y2)/2)1/2)不同,工况1d中的碎片床分布直径最大,为475.4 mm,如图14所示。在设计和优化堆芯捕集器时,应确定合理的烟囱顶盖垂直投影边长,使堆芯捕集器具备更好的性能。
图13 不同烟囱盖垂直投影边长下碎片床的y-z平面投影
图14 不同烟囱顶盖垂直投影边长下碎片床的最终分布直径
保持中央烟囱与周围烟囱间距为120 mm,烟囱顶盖垂直投影边长为90 min,设计烟囱顶盖倾角为15°~60°,形成第2组工况(工况2a~2g),如图15所示。通过数值模拟,得到碎片床的高度方差和高度归一化方差随堆芯捕集器烟囱顶盖倾角的变化结果,如图16所示。由图16可看出,随着烟囱顶盖倾角的增大,碎片床的高度方差和高度归一化方差呈先减小后增大的趋势。当烟囱顶盖倾角为45°时,堆芯捕集器托盘上形成的碎片床的高度方差和高度归一化方差达到最小值,分别为1 463.5 mm2和4.756×10-4。这表明,烟囱顶盖倾角为45°时,堆芯捕集器上的碎片床分布最均匀。需要指出的是,工况2b碎片床的统计方差小于工况2c,这是由于在工况2b中少量颗粒滞留在释放装置中,因此在统计方差时这些粒子被排除在外。
图15 工况2a~2g的堆芯捕集器烟囱顶盖倾斜角度示意图
图16 碎片床高度方差和高度归一化方差随堆芯捕集器烟囱顶盖倾角的变化
图17为下落颗粒在垂直于x-y平面并与x轴成45°夹角的投影。从图17可看出,除了工况2f和2g外,其他工况中碎片颗粒在烟囱顶盖上的滞留厚度均匀,这表明在下落过程中烟囱顶盖上的碎片颗粒数量处于动态平衡状态。在工况2a中,碎片颗粒几乎沿着烟囱顶盖向下滑落,没有观察到二次散射现象,且堆芯捕集器托盘上的碎片床分布直径也最小。在工况2e和2f中,可以观察到颗粒在堆芯捕集器烟囱顶盖上发生明显的二次散射。而在工况2f和2g中,烟囱顶盖上出现了明显的颗粒滞留现象,且二次散射的碎片颗粒数量也是该组模拟工况中最多的。
堆芯捕集器的烟囱顶盖对碎片颗粒的运动起到了减缓或反弹的作用,导致颗粒可能出现两种不同的后续行为。第一种运动行为是碎片颗粒被滞留在烟囱顶盖表面,并沿顶盖滑落,最终堆积在烟囱底部附近。另一种运动行为是碎片颗粒的二次散射,颗粒在烟囱顶盖表面发生反弹,继而落在更远的位置。在碎片床的形成过程中,一些颗粒滞留在烟囱顶盖表面形成堆积层,而另一些颗粒则被堆积层反弹形成散射层。当烟囱顶盖的倾斜角度增加到一定程度时,落到烟囱顶盖上的颗粒与从烟囱顶盖上二次弹射的颗粒发生大量碰撞,导致颗粒的散射效应逐渐被抵消。
需要指出的是,堆芯捕集器的最佳烟囱顶盖倾角与碎片颗粒的物性和机械性能相关。即在设计和优化具有相似结构的堆芯捕集器时,烟囱顶盖的倾斜角度应根据碎片颗粒的物性和机械性能来确定。通过优化堆芯捕集器烟囱顶盖倾角,可以使碎片床的形状和分布更加均匀,从而提升碎片床的热工水力性能。
保持烟囱顶盖倾角为15°,烟囱顶盖垂直投影边长为90 mm,设计中央烟囱与周围烟囱的距离为80~145 mm,形成第3组工况(工况3a~3l)。图18为该组工况碎片床颗粒的高度方差和高度归一化方差随堆芯捕集器烟囱间距的变化。结果显示,随着烟囱间距的增加,碎片床的高度方差和高度归一化方差呈整体下降趋势。当中央烟囱与周围烟囱间距为145 mm时,碎片床的方差达到最小值,表明工况3l中碎片床分布最为均匀。
图18 碎片床颗粒的高度方差和高度归一化方差随烟囱间距的变化
通过观察碎片颗粒的运动轨迹发现,中央烟囱与周围烟囱的间距对碎片颗粒的运动具有复杂的影响。当周围烟囱的位置恰好位于粒子二次散射的落点时,它们对碎片床的均匀分布起到积极作用。在该组工况中,碎片床的统计方差随着烟囱间距的增加呈整体下降趋势,表明所选取的烟囱间距并未符合周围烟囱位于颗粒二次落点的情况。然而,碎片床的统计方差在烟囱间距为95~115 mm的范围内出现了反复波动,由此推测粒子二次落点可能处于该范围内。但由于受计算工作量的限制,未对此猜想做进一步验证。
在设计和优化钠冷快堆的堆芯捕集器时,需要谨慎考虑周围烟囱的数量和位置。若周围烟囱处于碎片颗粒二次散射的落点,可改善碎片床的均匀分布;反之,周围烟囱的存在会阻碍碎片颗粒的分散。烟囱间距越大,周围烟囱的阻碍作用越小,这可通过图18中碎片床的统计方差得以证实。
本工作采用改进的无量纲系数离散元法,对钠冷快堆碎片床形成过程和堆芯捕集器的优化设计进行了数值研究。主要探究了堆芯捕集器的烟囱结构对碎片床的堆积形状和碎片分布的影响机制。模拟结果表明,碎片颗粒在堆芯捕集器烟囱顶盖上的二次散射明显改善了碎片床的形状和分布。在设计和优化类似结构的堆芯捕集器时,可以选择广泛的结构参数进行初步实验或模拟。根据本工作探究的碎片床形状和分布规律,缩小参数取值范围,并最终确定最佳设计参数。需要注意的是,堆芯捕集器的烟囱顶盖形状和烟囱高度也可能对碎片床的形状和分布产生影响,这方面的研究将在后续工作中展开。