具有外套管的释热元件表面温度求解

明 春，韩智杰，何晓军，初泉丽

(1.中国原子能科学研究院 反应堆工程技术研究部，北京 102413；2.国家核安保技术中心，北京 102401)

利用燃料组件导向管的内部空间开展材料或释热元件的辐照考验，是研究材料辐照性能及释热元件安全特性较常用的一种辐照试验方法。以典型压水堆燃料组件为例，燃料棒按照正方形栅格的形式排列，外侧为8根普通燃料棒，中心为导向管，内部插入释热元件棒。该结构中存在棒束通道外流流道-导向管-窄环缝内流流道-释热元件的复杂传热通道，其中，导向管内热流密度Q0未知，且内外流流道差异较大，实际计算时涉及到多个热边界的温度耦合问题，需联立求解多组传热微分方程，现有程序通常不会单独考虑这一特殊结构，难以计算出精确结果。工程上一般采用CFD软件进行仿真模拟，需进行三维建模，耗时较长、流程繁琐。因此，本文拟设计专用的耦合算法用于解决套管结构的传热问题，以便直接通过程序计算出符合燃料性能分析程序需要的套管结构温度数据。

实际计算时需要相应的实验数据进行传热关系式拟合。Weisman[1]和Markoczy[2]分别针对全堆芯和堆内有限区域给出了相应的传热关系式。对于环管流道，孙立成等[3]、曾和义等[4]和白博峰等[5]均开展过环管加热实验，得到了不同加热条件下环管传热的拟合关系式，Liu等[6]对环管传热的离心效应进行了深入研究。

本文拟基于以上文献，依据流场基本情况，通过迭代算法和热工水力基本理论，开发具有外套管的释热元件表面温度分布计算程序，并与燃料性能分析程序进行耦合，为具有外套管的释热元件的设计及安全评价提供技术支持。

1 几何模型

根据对称性，选择典型求解单元(图1虚线内范围)作为研究对象，不考虑组件所处具体堆芯位置对温度计算的影响。

导向管外侧流体所处流道为典型的棒束流道，本文将外流流道等效处理为圆形流道，则整个求解单元可视为1个多层套管结构，如图2所示，最外侧为燃料棒外壁，向内分别为外流流道、导向管、内流流道和释热元件。

为与性能分析程序计算节点相匹配，网格划分遵循1.5维的处理方式，首先将多层套管求解单元在轴向划分为若干个短圆柱段，对于每个轴向圆柱段，径向存在不同的传热部件。具体节点划分示于图3。

图1 求解单元结构示意图Fig.1 Structure of solving unit

图2 多层套管结构示意图Fig.2 Structure of multi-casing

图3 求解域节点划分Fig.3 Meshing of solution region

2 理论模型

2.1 冷却剂轴向传热

对于流体区域，已知两侧边界的热流密度，根据能量守恒关系有：

qΔl=cpqm(To-Ti)

(1)

其中：q为流体边界上线热流密度的代数和；Δl为换热边界长度，即轴向段长度；cp为比定压热容；qm为质量流量；To和Ti分别为轴向段出、入口温度，定性温度取轴向段流体平均温度。

2.2 外流流道对流换热

导向管外侧流体所处流道为典型的棒束流道[6]，在进行单向流分析时，需单独考虑通道形状对传热系数的影响。实际工程计算中，一般先采用Dittus-Boelter关系式[7]计算得到等效圆管的努塞尔数Nu∞,cir，然后进行系数修正得到棒束通道的Nu，即：

Nu=ψNu∞,cir

(2)

其中，ψ为修正系数。

对于正方形排列的棒束，考虑采用文献[8]给出的通用关系式计算ψ：

ψ=1+0.912Re-0.1Pr0.4(1-2.004 3e-B)

(3)

(4)

2.3 内流对流换热

估算可知，内流通道内Re<10 000，而Dittus-Boelter关系式仅适用于旺盛湍流换热，因此Nu需采用适用于过渡区的Gnielinski关系式[9]计算：

(5)

其中，下标f和w分别表示以流体温度和壁面温度作为定性温度。

根据白博峰等[5]的实验，环形管内的湍流换热相对于普通圆管内湍流换热有所加强。本文采用白博峰等进行的单侧加热环管实验得到的拟合关联式作为内流Nu的计算式：

Nu=0.023Re0.91Pr0.4

(6)

式(6)的适用范围为2 300

2.4 导向管导热

对于导向管区域，可视为无内热源的一维环形构件导热问题，导向管稳态传热方程[10]为：

(7)

其中：r为节点半径；T为节点温度；kc为导向管热导率。定性温度取导向管平均温度。

3 计算方法

稳态情况下，传热区域的边界上热流密度恒定，即导向管两侧热流密度Q0相同。通过求解Q0，可使内外流流道传热边界条件封闭，从而完成传热方程的求解。

本文所设计算法以轴向段作为基本求解单元，将Q0作为收敛标准，首先通过初始Q0以及传热相关方程逐步求解出外流流道、导向管和内流流道相应的节点温度，然后通过能量守恒反推出新的Q0，与初始值进行收敛判断，不收敛则重复上述过程直至收敛。Q0收敛后利用对流换热关系可求解出释热元件最外侧节点温度，完成当前轴向段温度计算，并将结果保存，作为下一轴向段计算的边值条件，最终完成全长温度计算。具体计算流程示于图4。

图4 轴向段温度求解流程Fig.4 Solution procedure of axial segment temperature

本文采用Fortran语言，依据此算法编译完成套管结构温度计算程序，用于对具有外套管的释热元件进行全长温度计算。

4 计算结果及对比

利用根据上述算法所建立的套管结构温度分布计算程序，采用典型压水堆燃料输入参数(表1)和ANSYS Fluent软件对程序进行对比分析。

表1 程序输入参数Table 1 Input parameters of program

4.1 计算结果

采用本文所建程序计算的套管结构内各区域的节点温度示于图5。由图5可见，外流流道与导向管外壁温度随高度的增加而上升，而导向管内壁、内流流道、释热元件外壁温度随高度的增加呈先上升后缓慢下降的趋势，温度峰值出现在出口附近。

图5 温度分布计算结果Fig.5 Calculation result of temperature distribution

4.2 结果对比

采用ANSYS Fluent软件对图1的3×3棒束结构建立三维模型，采用表1的边界参数进行热工水力仿真计算，并与本文程序计算结果进行对比，重点关注内、外流流道和释热元件包壳外侧的温度分布情况，以及导向管两侧的热流密度，结果示于图6、7。由图6、7可见，二者计算结果偏差较小，Fluent仿真模拟所得温度略高于程序计算值，相对误差在5%以内，本程序可在工程设计中提供一定参考。

图6 Q0本文程序计算值和Fluent模拟值的对比Fig.6 Comparison of Q0 calculated by code in this paper and simulated from Fluent

5 程序应用

套管结构温度计算程序改写为单独的计算模块后，可与传统性能分析程序相耦合，在性能分析计算开始前，调用该计算模块完成释热元件表面温度计算，结果作为边界条件在后续计算中使用。

图7 温度分布本文程序计算结果与Fluent模拟值的对比Fig.7 Comparison of temperature distribution calculated by code in this paper and simulated from Fluent

采用已耦合套管温度计算模块的燃料性能分析程序，进行完整的全寿期燃料性能分析，实现对燃料温度、变形、裂变气体释放等行为的研究[11-13]。通过假想算例，研究核电站全寿期内燃料性能的变化情况，尺寸参数及燃料棒平均线功率与表1一致，但释热元件线功率会随时间推移逐渐下降，因此其全寿期平均值设定为1 kW/m。

释热元件典型性能分析结果示于图8～10。由图8可见，寿期初芯块温度相对较高，温度分布曲线较陡峭，寿期末温度普遍降低，分布趋于平缓[14]。分析图9可知，前期由于裂变气体持续稳定释放，气体内压增大，而后期元件功率下降，裂变气体释放受到影响[15]，内压增长略有减缓。由图10可知，包壳与芯块的径向膨胀基本稳定，轴向伸长趋势则逐步减缓。

图8 寿期初与寿期末芯块轴向温度分布Fig.8 Axial temperature distribution of pellet during BOL and EOL

图9 寿期内元件内压的变化Fig.9 Variation of element inner pressure during lifetime

图10 寿期内元件力学变形的变化Fig.10 Variation of element mechanical deformation during lifetime

综上所述，基于本文算法修改的燃料性能分析程序在完整实现原有性能分析程序计算功能的前提下，可自行计算燃料棒外温度分布情况，并将传热方程求解所需边界条件传递给温度模块，保证后续运算流程正常进行，从而避免了对热工水力程序的依赖，提高了分析流程的独立性。

6 结论

本文针对套管传热结构设计了专门的迭代算法用于域内温度分布求解，并基于该算法开发了相应的Fortran程序模块，实现了具有外套管的释热元件包壳温度的求解。与大型CFD商用软件Fluent的仿真模拟结果比较，两者差异在5%以内，取得了较好的一致性。

耦合套管结构传热计算模块的燃料性能分析程序可自行完成温度边界条件的计算，独立完整地进行性能分析流程，减少了计算流程的复杂度。目前，该程序已实现工程应用，相关释热元件已入堆开展辐照试验。