基于多表面封闭系统网络法的钠冷快堆乏燃料组件内部传热数值分析

陈 翔，吴增辉，熊进标，程 旭，师 泰

(1.上海交通大学 核科学与工程学院，上海 200240；2.中国原子能科学研究院 反应堆工程技术研究部，北京 102413)

钠冷快堆乏燃料组件从卸出堆芯到进入乏燃料水池的转运过程中，会裸露在换热性能较差的气体环境中，自身的衰变热可能导致包壳过热失效，造成放射性泄漏事故，因而必须研究这一过程中乏燃料组件的传热特性，建立预测乏燃料组件最高温度的方法，特别是要对转运机构卡停事故的危险工况进行研究。

为获得不同工况下乏燃料组件的传热数据，国际上研究者进行了大量实验研究，如E-MAD[1]、SHATT[2]、REA[3]、Lovett[4]、SNL-LMFBR[5]等实验，其中SNL-LMFBR实验针对的是钠冷快堆乏燃料组件的厂外运输。除实验研究外，大量研究者利用FLUENT、StarCCM+等商用CFD软件对乏燃料组件换热特性进行了数值模拟研究[6-11]。

本文建立正六边形乏燃料组件的数值模型，开发数值分析程序研究乏燃料组件内部的传热特性，实现对组件内部每根燃料棒温度的预测。

1 传热模型

乏燃料组件内部传热有3种机制：导热、对流、辐射。本文模型仅考虑乏燃料组件径向方向上的传热，忽略导热和对流。

气体导热率较小而组件内部温度较高，在高温条件下辐射所占比重远大于导热所占比重，因此由导热带走的热量相比辐射可忽略。另外，因乏燃料组件内部流动阻力大，气体流动速度很小，导致对流传热带走的热量很少，故模型中不考虑对流换热。在忽略导热和对流的影响下，乏燃料组件仅靠辐射向外散热，因此本文模型相比实际情况更为保守。

本文采用的辐射模型基于多表面封闭系统网络法，以燃料棒外表面和组件盒内壁面构成封闭系统，并将封闭系统划分为若干计算节点。图1示出表面辐射示意图。每个节点入射到其表面上的辐射为G，其中一部分被吸收，占比α，称α为吸收比；另一部分被反射，占比ρ，称ρ为反射比。吸收比α和反射比ρ满足α+ρ=1。在漫灰表面的假设下，由基尔霍夫定律，吸收比α与表面发射率ε相等。

图1 表面辐射示意图Fig.1 Scheme of surface radiation

每个节点向外发射的辐射为J，称为有效辐射，包括自身的热辐射εEbA(Eb为黑体辐射力，A为表面面积)和向外反射的辐射ρG。因此每个节点满足方程：

Ji=εiEb,iAi+ρiGi

(1)

由玻尔兹曼定律，黑体辐射力Eb,i满足：

(2)

其中：σ为斯忒藩-玻尔兹曼常数；Ti为热力学温度。

角系数Xj→i是节点j发出的辐射入射到节点i的比例，其中，节点j发出的辐射为Jj，从节点j入射到节点i的辐射为Gj→i，于是，角系数为：

(3)

且角系数满足相对关系：

AiXi→j=AjXj→i

(4)

入射到节点i的辐射为各节点入射到节点i的辐射总和：

(5)

其中，N为节点总数。

结合式(1)～(3)、(5)，多表面封闭系统的辐射方程为：

(6)

每个节点满足能量平衡方程：

(7)

其中：cp，i为比定压热容；mi为质量；Qi为能量源项，由热源项Qsrc，i、辐射项Qrad，i组成：

Qi=Qsrc,i-Qrad,i

(8)

热源项由条件给定，辐射项为有效辐射Ji和入射辐射Gi的差值，满足：

(9)

式(6)、(7)构成传热模型的控制方程。

式(6)是有效辐射Ji的线性方程组：

(10)

程序中采用共轭梯度法求解。

图2 程序求解流程图Fig.2 Flow chart of solving procedure

2 实验验证

为验证数值模型的有效性，进行了37棒束正六边形乏燃料组件的模拟实验。图3示出组件横截面示意图。模拟组件内燃料棒呈三角形排布，长度为1 000 mm。37根燃料棒中30根为加热棒，内侧布置有加热丝，通过电加热模拟燃料棒的衰变功率。其余7根为测温棒，内侧不同高度处布置有热电偶，用以测量温度。7根测温棒处于图3中1、4的位置。

图3 组件横截面示意图Fig.3 Sectional view of assembly

为模拟组件建立数值模型，将燃料棒外表面和组件盒内壁面作为研究对象，划分成若干计算节点。每一根燃料棒外表面作为1个节点，这一类有37个节点；组件盒内壁面每一侧分成4个节点，节点分界点为外圈燃料棒中点在壁面上的投影点，这一类有24个节点，进一步地按照在组件盒内壁面一侧划分的长度不同分成两类：位于6个角上的12个节点，长度与棒间距、外圈燃料棒到内壁面距离有关；其余12个节点长度与棒间距相等。

上述节点布置考虑了结构的对称性，在此基础上将节点分为8组，如图3所示。前6组为燃料棒节点，第7组为边节点，第8组为角节点。特别地，第1组仅有1个节点，为中心棒节点。由于考虑了对称性，使得角系数容易计算得到[12-13]。简单的二维及三维几何结构，可得到解析公式，更复杂的几何结构需采用数值方法进行计算。模型中节点布置简单、对称，通过数学推导可得到任意两个节点间的角系数。

例如两根燃料棒间的角系数按间距从小到大排列，从节点1到节点2、3、5，角系数分别为：

(11)

(12)

(13)

其中，λ为棒间距与棒直径的比值。

其他角系数也有类似计算公式。每个节点所有角系数之和为1。数值分析程序中计算结果表明每个节点的角系数之和与1仅有计算精度的误差，间接证明了上述角系数计算公式的可靠性。

在实验工况的选择上，通过控制总加热功率，进行了功率为400、800和1 000 W 3组实验，分别测得组件盒内壁面温度为296、378和410 ℃。在数值分析程序中，给定内壁面温度和加热功率，计算得到每根燃料棒的温度。实验中7根测温棒没有加热功率，组件是非均匀加热，数值分析程序中设置成同样的非均匀加热条件。

除加热功率外，表面发射率在数值模型中作为另一组参数。表面发射率是材料的物理属性，然而由于技术原因难以测量，实验中缺乏相关数据。文献[14]表明，金属材料的表面发射率与表面粗糙度、氧化程度有关，而实验中由于高温条件会使得燃料棒表面氧化，因氧化程度不同导致表面发射率的不确定性。表1列出钢、不锈钢材料的表面发射率[14]。由表1可见，在抛光情况下表面发射率仅有0.1左右，而严重氧化后可达到0.8。以310型不锈钢为例，文献[15]表明表面发射率会随氧化程度、热处理温度的上升而从0.21上升至0.98。因此实验中表面发射率存在不确定性，数值模型中将表面发射率作为参数处理。在经过预计算后估计表面发射率在0.30～0.40之间，因此计算时选择0.30、0.35、0.40 3组工况，结合加热功率总计9种工况。

表1 钢材料的表面发射率Table 1 Emissivity of steel material

图4 数值模型计算结果与实验数据的比较Fig.4 Comparison between result of numerical model and experimental data

图4示出本文数值模型计算结果与实验数据的比较，纵坐标为节点温度T与壁面温度T壁面之差。由图4可见：表面发射率为0.40时的计算结果与实验值吻合较好，中心测温棒温度的偏差在10 ℃以内；外圈测温棒的计算结果均比实验数据要高，低功率条件下最为明显，偏差接近20 ℃。

表2列出表面发射率为0.40的工况下本文数值模型计算结果与实验数据的偏差。由表2可见：高功率时偏差很小，最高为6.3 ℃；低功率时偏差较高，特别是外圈测温棒，偏差达到18.3 ℃。这是由于低功率条件下组件整体温度较低，辐射传热所占比重减小，导热所占比重增大，而本文数值模型中忽略了导热的影响，导致此工况下计算结果过于保守，比实验数据偏大较多。而高功率条件下组件整体温度较高，辐射传热所占比重很高，满足本文模型的基本假设，因此模型在高功率条件下具有良好的可靠性。

表2 数值模型计算结果与实验数据的偏差Table 2 Deviation between result of numerical model and experimental data

3 Manteufel-Todreas双层传热模型

Manteufel等[16]建立了Manteufel-Todreas双层传热模型用以描述乏燃料组件内部的辐射传热和导热。模型在靠近边界处建立虚拟界面，将组件内部的封闭区域分为两个区域，称为核心区域和边缘区域，用不同公式描述两个区域的传热过程：

(14)

(15)

其中：S为形状因子，对于六角形排布的组件，取值12.836 5；Fpeak为轴向功率峰因子；La为加热段长度；Lc为截面周长；Fcond为导热率修正系数；Crad为辐射系数；Fcond，w为边界区域导热率修正系数；Crad，w，2为第二辐射系数；f为边界到内部的传热比；w为外圈棒中心与边界的距离；p′为列与列间的距离；kgas为气体导热率；Tm、Te分别为组件中心温度、虚拟界面温度；d为燃料棒直径。

按照实验中的几何尺寸和物理属性，式(14)、(15)的数值形式为：

Q=1.060 29(Tm-Te)+

(16)

Q=2.521 09×103(Te-Tw)+

(17)

用式(16)、(17)计算实验中的9种工况，与数值分析程序的结果进行对比。Manteufel-Todreas双层传热模型针对均匀加热的组件，因此数值分析程序中加热条件设置为均匀加热。

图5示出本文数值模型与Manteufel-Todreas双层传热模型计算得到的中心温度与壁面温度的差值，每种工况下数值分析程序的结果均比Manteufel-Todreas双层传热模型的结果偏高。本文数值模型中忽略导热的影响，而Manteufel-Todreas双层传热模型中考虑导热，因此本文数值模型相比更加保守，预测结果也比Manteufel-Todreas双层传热模型要高。

图5 数值模型与Manteufel-Todreas双层传热模型计算结果的对比Fig.5 Comparison between calculation result of numerical model and Manteufel-Todreas correlation

4 预测结果分析

实验方法和Manteufel-Todreas双层传热模型具有局限性：实验方法仅能得到测温棒处的温度；双层传热模型仅能得到中心处的温度。相比之下，本文数值模型能预测每根燃料棒的温度。

利用本文数值模型分析了不同工况下组件内温度的分布。图6示出非均匀、均匀加热条件下组件内温度的分布。由图6可见：非均匀加热条件下中心棒(节点1)并不是组件温度最高点，比临近节点温度要低，外圈测温棒(节点4)温度也比下一节点温度要低；而均匀加热条件下温度从内向外逐渐降低。依照式(8)热平衡时辐射项和热源项相等，非均匀加热条件下中心棒辐射项为0，那么中心棒温度需比其他棒温度要低以满足此条件，结果表现为中心棒温度比周围棒温度要低；同理外圈测温棒会比下一节点温度要低。

a——非均匀加热；b——均匀加热图6 不同加热条件下温度的分布Fig.6 Temperature distribution under different heating conditions

由图6还可见，最终功率分布的不同，导致温度分布有显著差异。以功率为1 000 W、表面发射率为0.40的工况进行比较，非均匀加热条件下，因测温棒(节点1、4)没有加热功率导致此处的温度相比均匀加热条件有显著下降，同时测温棒的热功率分配到加热棒处使得加热棒的功率上升，因此加热棒(节点2、3、5、6)温度有所上升。

以上预测结果表明，作为影响最高温度的重要因素，表面发射率、功率对组件温度的影响极其显著。表3列出非均匀加热条件下表面发射率为0.4时温度随功率的变化率。由表3可见，功率从400 W变化到800 W，功率每上升100 W温度上升10 ℃左右，而从800 W变化到1 000 W，功率每上升100 W温度仅上升7 ℃左右，变化率逐渐减小。

表面发射率下降会导致组件传热能力急剧下降。表4列出非均匀加热条件下功率为1 000 W时温度随表面发射率的变化率。由表4可见，温度随表面发射率的变化率非常大，且越靠近外圈变化率越大，而变化率会随表面发射率的增加而下降。

表3 非均匀加热条件下表面发射率为0.4时温度随功率的变化率Table 3 Temperature change ratio with power with emissivity of 0.4 under nonuniform heating condition

表4 非均匀加热条件下功率为1 000 W时温度随表面发射率的变化率Table 4 Temperature change ratio with emissivity with power of 1 000 W under nonuniform heating condition

5 结论

本文针对快堆乏燃料组件转运过程中处于气体环境时的传热过程建立了数值模型，基于多表面封闭系统网络法开发了数值分析程序，用于预测乏燃料组件内部每根燃料棒的温度。与模拟实验对比表明，高温条件下本文数值模型预测结果与实验结果偏差最高为6.3 ℃，证明了程序的可靠性。与Manteufel-Todreas双层传热模型相比，本文数值模型预测结果整体高于Manteufel-Todreas双层传热模型预测结果，表明本文数值模型更具保守性。在此基础上分析了加热功率、表面发射率、均匀加热、非均匀加热对温度分布的影响。结果表明，温度随功率的上升而上升，表面发射率下降会导致组件传热能力急剧下降。同时，均匀加热时燃料棒温度从中心依次向外降低；非均匀加热时燃料棒温度分布会发生偏移，最高温度并不位于中心位置。