基于COMSOL的板状燃料组件流道堵塞事故研究

袁显宝 刘 超 谭 伟 张永红 张彬航 毛璋亮 周建军 唐海波 刘芙蓉

(1.三峡大学 机械与动力学院，湖北 宜昌 443002；2.湖北省水电机械设备设计与维护重点实验室，湖北 宜昌 443002)

堆芯流道堵塞是严重的堆芯事故之一.由于板状燃料组件的设计紧凑，当发生燃料肿胀、堆内材料碎片或者异物进入堆芯循环等异常事件时，就可能引发堆芯堵流事故.此时，反应堆内冷却剂的流动受到阻碍导致传热恶化，燃料释放热量的滞留将引起组件温度升高并威胁燃料包壳的完整性，造成安全事故，导致反应堆停止运行[1].

对于窄缝结构内部温度场和流场的研究，由于测量技术限制，实验实施困难，相关研究主要采取数值模拟的方法[2].计算流体动力学(computational fluicl dynamics, CFD)方法可构建精细的三维几何结构，直接求解研究对象的能量守恒方程；并对堵塞区域进行合适的模拟，得到清晰的可视化结果，在反应堆热工水力研究中得到了广泛应用.宋磊[3]等使用Fluent软件对板状燃料组件进行了三维CFD计算，求解了单流道95%堵塞和全部堵塞造成的燃料组件内的热工水力变化.董化平、樊文远[4-5]等利用Fluent研究了板型燃料组件和多层环形板状燃料组件在发生入口堵塞后物理场的变化情况，探索了使用动网格技术模拟堵塞面的方法在研究该类事故中的应用效果.Amgad Salama[6]对典型材料测试堆板状燃料组件发生弯曲堵流现象进行模拟，研究堵塞程度导致的堆型流场变化和温度变化.

目前国内对板状燃料组件发生堵塞的研究相对较少，本文基于COMSOL Multiphysics有限元软件的热流耦合技术，开展了板状燃料组件单流道堵塞工况下的流体流动和流固传热数值研究，对板状燃料组件流道堵塞事故的预防和事故严重性评估具有一定的参考价值.

1 数值模型

计算模型基于国际原子能机构IAEA(International Atomic Energy Agency)10 MW 轻水冷却和慢化的理想化池式材料测试堆(MTR，Material Test Reactor)，堆芯参数[3、6]见表1.

表1 MTR堆芯参数

1.1 几何描述

MTR堆芯共包含21盒组件，其中1个标准燃料组件包含23块燃料板，燃料板以平行等间距的方式排列.组件俯视图如图1所示.燃料板由燃料芯体和

图1 燃料元件结构(单位：mm)

外部覆盖的包壳组成，冷却剂沿着Z轴负方向通过板间缝隙并带走燃料芯体所释放的热量.模型的上方和下方各设置高度为300 mm的空腔，使流体流动更符合堆芯运行时燃料组件周围的冷却环境[3-5].表2给出了模型几何参数.

表2 燃料元件的几何参数

文献调研结果表明，板状燃料堵塞事故主要影响堵塞流道及相邻的两个流道[7]，故相关研究常选取组件边缘起2到3块燃料板及相邻冷却剂流道作为研究对象[1].考虑到板状燃料布置和几何的对称性，选取了单组件边缘2块板和3个通道的半部分进行建模.

当组件内流道发生堵塞时，由于阻塞物体的形状和厚度是不可预测的，根据典型堵流事故实例[8]，堵塞原因以堆外异物掉入为主，极有可能将整个流道入口完全覆盖.依照樊文远[4]提出的4种假定的入口堵塞形式，采取在流道入口设置刚性无厚度薄面的方式来模拟阻塞工况，属于流道堵塞工况中的中心堵塞.图2中左侧图形为模型的原本几何形状，考虑到燃料组件纵横比较大，视觉上缩放为图2右侧容易辨识的三维几何图形.以下腔室一顶点为原点构建三维模型，其俯视图的a、b、c、d4个顶点与图1平面图中所选的研究区域相对应.

图2 组件流道发生堵塞的几何模型

1.2 数学模型

假设冷却剂为牛顿流体，在流道中的流动为无相变湍流流动.根据Daxin Gong[9]对不同湍流模型在板状燃料组件计算适用性的研究，选取realizablek-ε两方程模型作为求解湍流流场的数学模型.该模型是k-ε两方程模型的修正形式，能对范围较广的湍流流动现象进行可信模拟，同时具有良好的鲁棒性和计算经济性，通过引用两个新的变量湍流动能k与湍流耗散率ε来求解湍流粘度μt.使用壁函数来近似替代近壁区域流体的法向速度变化.对于湍流区域，控制湍流动能k的输运方程和湍流动能耗散率ε的输运方程组为：

湍流动能方程：

Gk+Gb-ρε-YM-Sk

(1)

扩散方程：

(2)

其中：

(3)

式中：ρ为流体密度；μ为流体动力粘度；ν为流体的运动粘性系数；uj为流体在xj方向的流动速度；Gk为平均速度梯度面产生的湍动能；Gb为受浮升力而产生的湍动能；YM为可压缩湍流中脉动膨胀对整体湍流扩散率的贡献，在不可压缩流体的计算中通常忽略掉；C1ε为常量，计算中取1.44；C2为常量，计算中取1.9；C3ε为浮力对ε的影响系数；σk为k的普朗特数，计算中取1.0；σε为ε的湍流普朗特数，计算中取 1.2；Sk为源项；Sε为源项.

1.3 网格划分

由于板状燃料间隙狭窄，仅有2.23 mm，燃料板厚度仅有1.27 mm，而模型在高度方向整体长1 200 mm，纵横比较大，对网格精度要求较高.冷却剂在由上腔室进入流道和从流道中流出时，湍流强度高，速度变化剧烈，因此需要对相应区域的网格加强细化.

参考Amgad Salama[6]对燃料板组件堵流事故进行研究时二维计算模型的网格划分方式，采取结构化网格和非结构化网格相结合的方式对模型进行网格划分，如图3所示.对燃料板及流道内冷却剂区域进行划分时采取结构化网格，上下腔的流体区域采取非结构化网格，边界层设为3层.进行无关性验证后，网格包含400万个计算单元，整体平均质量为0.896 6.

图3 板状燃料元件模型的网格划分

1.4 边界条件

设置湍流模块和流体与固体传热模块，湍流模块中依次设置入口边界条件为速度入口、出口边界条件为0、压力出口以及壁条件为壁函数.为了简化模型设置了两个对称边界.采用内壁条件模拟堵塞面，该边界位于Ch1通道顶端，改变该边界的面积占比来控制堵塞程度.对称面、堵塞面位置详见图2.流体与固体传热模块中分别选定燃料板为固体域、冷却剂所在区域为流体域，流入温度为311.15 K，出口边界和对称边界与湍流模块相一致.设置燃料芯体为热源，释热功率按体积均匀分布.建模所需具体参数和初始值设定见表3.

表3 边界条件和初始值[3]

2 结果分析

2.1 正常运行状态

正常运行时，冷却剂自顶向下流动，在燃料板顶面受到阻碍形成3条分流，分别流经3个狭窄通道(Ch1、Ch2、Ch3)，在下腔室汇聚，最终从下腔室底面流出.图4中云图为模型中分别位于上腔室、燃料区、下腔室3个区域XY截面上冷却剂的流速分布.

图4 正常运行时冷却剂在不同截面的流速分布

图5为对应图4中截线Y=16.625 mm上冷却剂流速在X轴的分布.冷却剂以4.1 m/s速度进入上腔室，即将进入流道时，流速分布较平坦，右侧的近壁区流速出现衰减；冷却剂受燃料板阻隔后分流，在流道内流速分布趋于一致，流速峰值在流道中心，达到7 m/s；经流道流出后，右侧近壁区的流速在流动方向上衰减缓慢，相对流速较高，该结果与参考文献[7]一致.

图5 组件不同高度截线上冷却剂流速随X轴变化

图6所示云图为XZ面上的温度分布，位置为Y=16.625 mm.由于各流道冷却剂流量大致相同，在各燃料板功率分布一致的条件下，组件内部的温度场分布基本对称.燃料板主要以对流换热方式向冷却剂转移自身产热.在流动方向上，随着冷却剂被加热，与燃料板温差减小，换热效率逐渐降低，冷却作用逐渐减小，燃料板的温度逐渐上升，接近底部时达到342 K；在X轴方向上，燃料板与冷却剂温差较大，包壳附近温度梯度较大，但流道中心大部分区域温度较低，冷却剂对燃料板释热仍有一定裕量.

图6 Y=16.625 mm面的温度分布

2.2 堵塞事故状态

堵塞会对冷却剂的正常流动造成较大影响，引起流量重新分配.如图7(a)所示，在流道入口附近，堵塞面的存在迫使原本将流入Ch1的冷却剂向旁边流道转移.流道出口附近，堵塞流道Ch1到下腔室的过渡区域发生较为明显的漩涡和回流现象，如图7(b)所示.经过对比各流道入口横截面上的平均流速，得到表4关于不同运行工况下各冷却剂流量占比.

图7 堵塞引起流量转移和回流

表4 3种工况下各流道流量占比 (单位：%)

堵塞发生后，流量再分配和流道内漩涡的存在使得相应区域的温度发生明显变化，如图8所示.由图可知，未发生堵塞时，两燃料板温差较小，Ch1内冷却剂平均温度与入口温度相比无明显升高；发生95%堵塞后，Ch1流道的冷却剂流动速度减慢，燃料板热量转移和导出效率降低，热量滞留在Ch1流道内导致冷却剂温升明显.左侧燃料板出现不连续的高温区域，Ch2、Ch3流道及右侧燃料板受影响较小；全堵塞时，左侧燃料板温度出现了连续的高温区域，Ch1流道内冷却剂温度大幅升高且与左侧燃料板温度十分接近，但右侧燃料和Ch2、Ch3流道的温度变化仍不明显.

图8 组件中轴线的温度分布

图9为不同工况下三维截线(面Y=16.625 mm与面Z=600 mm的交线)上的温度分布.在正常工况下，流经3个流道的冷却剂流量基本相同，燃料板向两侧传递的热量基本相同，因此温度基本呈对称分布，Ch3流道由于右侧无热源，整体温度稍低于Ch1.95%堵塞工况下，Ch1流道的流量急剧下降，整个流道内的温度大幅升高，部分区域达到了正常运行时燃料板处温度340 K.堵塞引起的流量降低和回流现象使温度峰向左侧小幅度偏移，但影响范围仅限Ch1流道和左侧燃料板，此时Ch2流道的温度上升现象并不明显.由于堵塞使流经Ch1的流量向Ch2和Ch3转移，使得右侧燃料板与冷却剂之间的换热更加充分，温度小幅降低.

图9 组件中轴线的温度分布

全部堵塞工况下，Ch1流道的冷却剂几乎无法提供冷却作用，燃料板温度进一步升高，达到了360 K.但同95%堵塞状况相似，全堵塞工况对右侧燃料板传热无明显影响.中轴线上的温度分布状况与参考文献[3]中Fluent计算所得结果(E)偏差较小.

表5给出了3种工况下，组件不同区域的最高温度.全堵塞工况相对于95%堵塞工况的各部分区域最高温度均有上升，包壳的最高温度为362.74 K，冷却剂最高温度360.45 K，均未超过环境压力0.17 MPa时的饱和温度388 K[7]，维持在安全范围内，因此不会在燃料板包壳表面产生过热沸腾现象.

表5 堵塞前后不同区域中的最高温度 (单位：K)

3 结 论

以IAEA 10 MW MTR堆为对象，采用COMSOL程序对其标准燃料组件部分结构进行三维建模并求解了堆芯发生流道堵塞时热工水力参数变化情况，得出以下结论：

1)稳态运行时，组件在X方向温度呈周期性对称分布，沿冷却剂的流动方向温度逐渐升高，接近底端时到达温度峰值，约为342.73 K.

2)发生95%堵塞时，受阻塞影响冷却剂流量将重新分配，堵塞流道Ch1内的冷却剂将形成漩涡并影响传热.由于冷却剂仍保持流动，传热恶化的程度有限.完全堵塞时，左侧料板单侧丧失冷却，温度迅速升高，燃料包壳的最高温度达到了362.74 K，但未超过冷却剂饱和温度，不会引起沸腾.单流道堵塞影响区域仅限于堵塞流道及相邻的燃料板，相邻流道温度升高不明显.

3)完全堵塞时，堵塞流道内的冷却剂温度可达360.45 K，考虑到每块燃料板仅有两侧流道，流道间隙较狭窄，存在相邻两流道同时受到阻塞的可能性.若发生相邻两流道均被堵塞的情况，整块燃料板释放的热量将无法顺利导出，温度升高程度更严重，可能对组件结构完整性产生较大威胁.