钠冷快堆六角形组件换热特性分析

2021-02-03张东辉刘一哲

原子能科学技术 2021年2期

师泰，张东辉，刘一哲

(中国原子能科学研究院，北京 102413)

示范快堆是我国自主设计建造的一座池式钠冷快中子反应堆，反应堆采用六角形组件作为燃料，燃料在反应堆运行过程中浸泡在导热性能较好的液态金属钠中，在换料过程中会通过换热性能较差的氩气空间。因此，保证燃料组件在氩气环境中安全性是反应堆安全设计的重要组成部分[1]。快堆燃料组件初始释热较大，一般采用乏燃料组件装入堆内储存阱中放置1～2个换料周期的方式降低其衰变热功率。换料开始后，组件由换料系统提升到反应堆外部的转运室中，转运室为氩气环境，长时间悬停在转运室可能导致组件超温，甚至破损，组件破损熔化会导致放射性物质进入环境[2-3]。因此，研究燃料组件在氩气环境下的换热特性，保证反应堆换料工程中组件的安全极其重要。

目前国内外针对六角形组件在气体环境下的换热特性研究甚少，美国桑迪亚国家实验室[4]开展了模拟水平运输过程快堆217棒乏燃料组件的换热特性试验，进行了低加热功率范围在不同气体介质中的传热试验，得到了组件内部温度分布并拟合了热导率经验关系式[5]。Alyokhina等[4-6]建立含多个快堆乏燃料组件的储存桶的三维模型，研究环境温度变化对贮存桶内温度分布的影响，结果表明外部环境温度变化对贮存桶内温度分布无明显影响，同时解决了反向共轭传热问题。

本文采用三维CFD软件研究37棒乏燃料组件在充满氩气环境的转运室悬停时，燃料组件依靠自然循环方式冷却的稳态温度场，同时采用37棒模拟组件实验进行验证。

1 实验简介

该实验模拟37棒组件在氩气环境下的换热，试验装置如图1所示。试验段悬挂在压力容器中，通过对试验段组件进行电加热模拟组件的发热，组件内布置有热电偶，通过热电偶测量组件内的温度分布，从而得到组件在不同功率下的局部温度分布。压力容器采用三段式设计，氩气容器采用真空泵抽真空后注入氩气，氩气压力稳定在(0.10±0.01) MPa，为满足容器壁面稳定的温度边界条件，在压力容器外壁焊接有控制边界温度冷却管路。压力容器内径为500 mm，分3段设计，下段高度为400 mm，中段高度为1 200 mm，上段高度为800 mm。

图1 试验装置示意图Fig.1 Schematic of experimental apparatus

图2 试验段Fig.2 Test section

37棒模拟组件试验段竖直安装于压力容器中段中心位置，试验段由棒束和组件盒组成，如图2所示。棒束呈六角形排列，由30根加热棒和7根测温棒组成。加热棒的内部结构如图3所示，棒内使用导热性能较好的MgO粉末填充压实。根据导电率随温度的变化关系，将处于对称位置的加热棒组合为5组，即H1～H5，如图2b所示。每组加热棒并联后与程控电源连接，独立调节电源电压，保证加热功率均匀。

组件盒为正六边形结构，内对边距为44 mm，厚度为3 mm。元件棒直径为6 mm，棒中心距为6.95 mm，总长度为1.2 m。

图3 加热棒的结构Fig.3 Structure of heated rod

中心测温棒(M1)内安装4只热电偶，外圈测温棒(M2～M7)分别安装2只热电偶，测温棒不具备加热功能，测温棒的测点位置如图4所示。

图4 M1～M7内热电偶位置Fig.4 Thermocouple locations in M1-M7

2 计算程序及模型

2.1 计算程序

本研究采用 CD-asapco公司开发的STAR-CCM+软件进行计算，该软件采用最先进的连续介质力学数值技术开发的新一代CFD求解器。它搭载了CD-asapco独创的最新网格生成技术，可完成复杂形状数据的输入。它的主要特点在于条件设置与后处理方面，如同时打开两个算例，某些条件设置仅需复制粘贴即可互相转换，省去了再输入1遍的时间，尤其是在气体吸收系数各波段的波长范围录入以及吸收系数录入等数据量大的设置中，为用户节省了大量时间。另一方面，STAR-CCM+在多面体网格生成方面有较大的优势，尤其是应对复杂几何结构的计算，如燃料棒束、堆芯等复杂几何体[7-9]。

2.2 关键模型分析

1) 自然循环能力估算

热能传递的方式有3种：热传导、热对流和热辐射。乏燃料在氩气环境下没有强制对流，组件内部靠自然循环建立流动特性，由于组件内部间隙较小，自然循环能力低，可忽略不计。自然循环能力估算如下。

假设组件与氩气充分换热，氩气出口温度与组件温度一致，元件棒表面温度为700 ℃，氩气在组件内的自然循环流量计算公式[10-11]为：

(1)

其中：qm为流体质量流量，kg/s；ρ为密度，kg/m3；g为重力加速度，m/s2；A为流体的横截面积，m2；αv为体膨胀系数，K-1；tw为组件表面温度，℃；ta为氩气环境温度(与环境温度一致，50 ℃)；δ为元件棒之间的间隙；ν为气体的运动黏度，m2/s。

氩气自然循环所带走的热量为：

Q=Δt·qm·cp

(2)

其中：Q为热量，W；Δt为氩气温度变化量，℃；cp为比定压热容，J/(kg·℃)。计算得到Q为17.187 W，相对于400～1 000 W组件功率可忽略不计。

2) 绕丝的模型简化

由于组件内部流动特性较小，绕丝对组件换热特性的影响主要由于绕丝较好的导热性能，而绕丝对组件内部流动的影响可忽略。因此，CFD建模中采用无绕丝的模型，绕丝对组件内部导热的影响采用等效导热的方法分析，无绕丝的等效导热方法相对于真实组件的计算结果保守。采用等效热阻方法得到等效导热系数[12-14]如下：

(3)

其中：λt为等效导热系数；λa和λs分别为氩气和钢的导热系数,W/(m·K)；Aa和As分别为氩气和钢的等效面积,m2。氩气的导热系数为0.036 W/(m·K)，计算得到等效导热系数为0.051 W/(m·K)。

3) 辐射模型

S2S(surface-to-surface)辐射模型适用于计算封闭系统中各反射面之间的辐射传热，并且模型主要参数运行在计算迭代步数开始前预处理完成，因此计算具有较高的效率，故本文选择S2S模型[15]。

对空间两个微元面而言，其辐射换热如图5所示，由dS1面发射被dS2面吸收的辐射总功率为：

P1-2=i′1dS1cosβ1(dS2cosβ2/L2)

(4)

其中：P为功率；i′为辐射强度，W/(m2·sr)；S为面积；L为特征长度；β为空间角度。

图5 微元面辐射换热Fig.5 Micro-element radiative heat transfer

两微元面之间的辐射角系数可由下式计算获得：

(5)

其中，F为角系数。当两个黑体表面在相同温度场中达到热平衡时，其相应的辐射角系数应满足如下关系式：

dFi-jdSi=dFj-idSj

(6)

据此拓展到空间有限大两平面之间时，有：

(7)

同理可知：

(8)

2.3 边界条件设置及计算工况

乏燃料组件在转运室的换热是一个自然对流冷却的过程，由于自然循环能力较弱，因此忽略对流的影响，组件主要通过导热和辐射散热，计算为稳态计算。边界条件主要包括热源边界、固体边界和外壁面边界。

1) 热源边界，乏燃料组件在出堆后组件功率轴向分布可近似为恒定值，因此元件棒功率采用体积功率的输入参数，对30根元件棒轴向0.2～1 m段赋予相同的体积功率。实验中加热棒为30根元件棒，因此400、800和1 000 W功率下体功率分别为5.89×105、1.18×106、1.47×106W/m3。

2) 固体边界，乏燃料组件内部元件棒与氩气和氩气与组件盒之间为接触界面边界条件。固体导热系数依据指定的导热系数进行计算，固体之间无接触热阻。辐射传输采用氩气内部辐射传输。

3) 外壁面边界，外壁面边界条件主要包括组件盒外壁和上下端面。该边界条件采用壁面边界条件，由于实验过程中组件出现部分氧化现象，组件盒外壁的发射率如图6所示，表面对流换热系数均采用5 W/(m2·K)，环境温度为常数50 ℃。上下端面对组件换热影响较小，故采用绝热壁面边界条件。

图6 归一化发射率分布Fig.6 Normalized emissivity distribution

CFD计算与实验的工况相同，计算输入参数为：氩气导热系数，0.051 W/(m·K)；组件盒和元件棒导热系数，16 W/(m·K)；表面换热系数，5 W/(m2·K)；环境温度，50 ℃；功率，400、800、1 000 W。气体的等效导热率为0.051 W/(m·K)。

2.4 网格划分及收敛性

多面体网格相对于四面体网格能在复杂几何体及边界层、长通道或小间隙的区域可在较少的控制体的基础上达到较高的精度。同时，由于多面体网格有很多邻居单元，所以能精确计算控制体的梯度。但多面体网格可能造成较大的计算量和内存需求。由于组件尺寸跨度大、结构复杂，采用多面体网格可提高计算精度。整体网格如图7所示，网格质量列于表1。表1中，网格在拓扑上有效，并且没有负体积，总网格数为5 035 336。

图7 37棒组件网格图Fig.7 37 rods assembly mesh

表1 网格质量Table 1 Mesh quality

本文分别采用200万、400万、500万、900万和2 200万网格进行了测试计算，中心棒600 mm位置温度如图8所示。结果表明，500万网格与2 200万网格对应的组件中心棒600 mm位置温度的相对误差为0.29%，能满足要求，故采用500万网格进行分析。

由于乏燃料组件仅考虑导热和辐射，组件盒内氩气采用固体模型，计算收敛性主要参考能量残差和局部温度恒定。图9示出计算迭代过程中的残差和中心棒800 mm处的温度随迭代步数的变化，图中可看到迭代步数4 500步之后组件温度基本不变，能量残差小于10-5，故判断计算收敛。

图8 网格数与测量值的关系Fig.8 Relationship between number of grid and measured value

图9 能量残差和温度随迭代步数的变化Fig.9 Variation of energy residual and temperature with number of iteration step

3 实验和计算结果分析

乏燃料组件温度直接关系到乏燃料组件操作的设计及安全性，尤其是燃料元件包壳的最高温度，是关系到组件能否包容放射性的重要因素。数值模拟结果与实验测量值相对误差较小，能满足CFD模拟实验的精度要求。由于采用的填充材料为氧化镁，该材料导热性能较好，实验中，中心棒600～800 mm的温差较小(<6 ℃)，可近似认为实验和计算的中心棒600、650、800 mm 3个位置的最高温度为组件最高温度。

实验和数值模拟结果表明中心测温棒的温度最高，且最高温度出现在600～800 mm之间，400 W功率下最高温度为440 ℃，800 W功率下最高温度为581 ℃，1 000 W功率下最高温度为631 ℃。实验和数值模拟不同功率下中心棒M1最高温度对比如图10a所示，实验和数值模拟不同功率下盒壁最高温度对比如图10b所示，实验和数值模拟中心棒M1的温度分布对比如图11所示。

图10 不同功率下中心棒M1和盒壁最高温度对比Fig.10 Comparison of the highest temperatures for central rod M1 and component box at different powers

图11 不同功率下中心棒M1温度对比Fig.11 Comparison of temperature for central rod M1 at different powers

不同功率下中心棒M1温度CFD计算值和实验值相对误差对比如图12所示。由图12可知，实验数据和CFD模拟的相对误差较小，最大的相对误差为6%，出现在距离上下端较近位置，这是由于在实验中采用钢条对实验装置进行了固定，钢条较好的导热性能对实验产生较大的影响。CFD计算中最高温度能较为精确地反映实验中最高温度，中心棒600 mm和650 mm位置的温度相对误差不超过1%。

图12 不同功率下中心棒M1温度的相对误差Fig.12 Relative error of temperature for central rod M1 at different powers

实验和数值模拟结果表明外圈组件温度较中心棒组件温度低，由于实验和CFD计算中外圈测温棒不同位置的测温点不在同一元件棒上，由于不同元件棒相对于中心棒的角系数有一定的区别，辐射换热传递的热量不同导致温度出现较小的波动。外圈组件的最高温度出现在600～800 mm之间，CFD模拟中400 W功率下最高温度为407 ℃，800 W功率下最高温度为537 ℃，1 000 W功率下最高温度为587 ℃。实验和数值模拟外圈测温棒M2的温度分布对比如图13所示。不同功率下外围组件温度CFD计算值和实验值相对误差对比如图14所示。由图14可知，实验数据和CFD模拟的相对误差较小，最大的相对误差为3%。

图13 不同功率下外圈测温棒M2温度对比Fig.13 Comparison of temperature for outer test rod M2 at different powers

图14 不同功率下外围测温棒M2温度的相对误差Fig.14 Relative error of temperature for outer test rod M2 at different powers

实验和数值模拟结果表明盒壁温度较燃料元件温度低，由于实验和CFD中盒壁测温点处于六角形各边的中心，测量温度为盒壁的最高温度，盒壁温度在6个角位置较低。盒壁最高温度出现在600～800 mm之间，CFD模拟中400 W功率下最高温度为285 ℃，800 W功率下最高温度为373 ℃，1 000 W功率下最高温度为405 ℃。实验和CFD模拟盒壁温度分布对比如图15所示。不同功率下外围组件温度CFD计算值和实验值相对误差对比如图16所示。由图16可知，实验数据和CFD模拟的相对误差较小，最大的相对误差为9%。