Study on the Cross-flow Characteristics of Wire-wrapped Fuel Bundle Channel in CiADS

铅铋合金(Lead-Bismuth Eutectic, LBE)有很多优良特性，如中子经济性良好、导热系数高好、沸点高，且不与水发生剧烈化学反应等。由此，加速器驱动嬗变研究装置(China initiative Accelerator Driven System, CiADS)^[1-3]采用LBE为次临界反应堆的堆芯冷却剂。同时，CiADS拟采用绕丝作为燃料棒的定位部件，可增强冷却剂在燃料组件内的横向交混，展平燃料组件内冷却剂的温度分布，进而提高反应堆的安全性和经济性。因此，研究绕丝引起的横流交混特性对CiADS燃料组件的热工水力学分析具有重要意义。

对带绕丝燃料组件流动特性的研究多基于计算流体动力学方法。由于绕丝几何模型复杂，在CFD方法中研究者通常采用非结构化的网格研究横流交混特性。Ahmad等^[4]使用雷诺时间平均模型(Reynolds Averaged Navier Stokes, RANS)分析了绕丝对冷却剂流动的影响，发现绕丝能使冷却剂的交混更加充分。Merzari等^[5]对7棒束绕丝燃料组件开展了大涡模拟(Large Eddy Simulation, LES)和RANS研究，发现k-$\omega $ SST湍流模型的计算结果与LES的结果最接近。Song等^[6]将CFD模拟结果与粒子图像测速法(Particle Image Velocimetry, PIV)的实验数据进行了对比，为CFD提供了横向速度基准验证数据。Li等^[7]研究了流道结构对流量分配的影响，发现采用更小外圈内径的绕丝和减小对边距可以对不同子通道中流量进行重新分配，能使出口温度分布更加均匀。Fan等^[8]也基于非结构化网格对CiADS燃料组件进行CFD模拟，发现子通道间隙横向速度特性的周期性。但是，非结构化网格会损失小尺度的涡流细节，且存在网格数量多和部分网格质量低的问题。而结构化网格质量更高，可获得更精细的流动细节，因此有必要采取该方法对CiADS带绕丝燃料棒束横流交混特性进行研究。

本论文将采用结构化网格和CFD方法对CiADS燃料棒束通道横流特性进行研究，分析结构化网格与非结构化网格对流动特性的影响，将基于更为准确的计算结果对横流特性进行定量分析 。

CiADS堆芯包含30盒燃料组件，每盒燃料组件由呈三角形排布的61根棒组成，其中60根为燃料棒，剩下1根为位于中心的不锈钢棒^[1]。相邻燃料棒之间利用螺旋绕丝进行定位，绕丝沿流动方向呈右手螺旋分布。燃料棒外径为13.1 mm，绕丝直径为2.0 mm，相邻燃料棒间距为15.1 mm，绕丝螺距为144 mm，对边距为122.5 mm。本研究使用周期性边界条件进行计算，简化为单螺距流道^[9]，有利于减小计算量和提高网格密度，以获得更精确的结果。单螺距的燃料棒束流道三维模型如图1所示。

划分计算网格是CFD计算中较为重要且费时的一步，网格的划分方式、密度和网格质量对计算精度有很大影响。划分带绕丝燃料棒束通道的网格，一般需要进行简化处理：即将绕丝往燃料棒方向缩进微小距离。同时，非结构化网格一般还需要在绕丝和燃料棒间的尖角进行倒圆角处理以避免线接触，影响网格质量。Hamman等^[10]发现倒圆角简化处理虽然对模拟结果整体影响不大，但牺牲了局部流场细节。

本研究使用ANSYS ICEM^[11]软件进行结构化网格划分，将绕丝沿燃料棒缩进0.1 mm，无倒圆角处理。结构化网格如图2(a)所示，分成三个步骤进行划分：(1) 划分燃料棒和绕丝的包裹流道；(2) 划分间隙流道；(3) 进行网格合并。其中：燃料棒和绕丝的包裹流道采用二维面网格拉伸的方法，以燃料棒中心线为轴进行旋转拉伸，当拉伸长度到达144 mm时，面网格恰好旋转一周。间隙流道采用构建三维block进行线和面关联的方法生成网格。轴向方向上，燃料棒周围流道与燃料棒外部流道的交界面处设置多条螺旋线作为关联线；横截面上只需要取进出口两面进行线面关联。ICEM中网格质量检测显示，99%以上网格的质量在0.32以上，网格质量较低的地方为流动死区，整体网格质量满足计算要求。将两部分网格合并、导入CFD求解器软件ANSYS Fluent^[12]后，网格重合面设置为Interface边界条件以进行边界面关联。边界层第一层厚度通过估算y+确定，即使无量纲壁面距离满足RANS湍流模型的应用条件。本研究中雷诺数为20 000工况下，壁面最大y+为2，在湍流模型的适用范围以内，边界层厚度设置合理。

此外，为研究结构化网格与非结构化网格之间的差异。本研究使用STAR-CCM+^[13]软件进行了多面体网格划分。多面体网格对几何模型适应性强、计算效率高，计算结果与四面体网格非常接近，可以达到较高的精度^[14]。据此本研究将其选用为非结构化网格代表。网格物理模型采用的简化方式为缩进距离0.1 mm，绕丝与燃料棒之间缝隙倒圆角半径0.25 mm。网格生成的主要参数为：基础尺寸为0.6 mm，边界层4层，网格体积增长率为1.2，雷诺数为20 000工况下，壁面最大y+小于2。网格示意图如图3所示。

数值模拟方法的精度由高到低依次为：直接数值模拟(Direct Numerical Simulation, DNS)、LES和RANS；然而，DNS和LES方法对网格质量和网格密度的要求更高，计算量大，在工程上实用性不高。Dovizio等^[9]的研究表明：RANS的计算结果与DNS、LES的计算结果吻合良好，可使用RANS方法代替DNS、LES方法进行研究，节约计算成本。Merzari等^[5]的研究表明：RANS模型中，k-$\omega $ SST模型^[15]计算结果与实验对比吻合最好，这是因为k-$\omega $ SST模型在保持了k-$\omega $模型对剪切流强的捕捉能力的同时，也能对充分发展流进行较好的预测。因此，本研究选用k-$\omega $ SST模型进行迭代计算。

本研究模拟了雷诺数从5 000到20 000之间不同工况的流动特性，选用周期性质量流量作为边界条件，采用SIMPLE^[12]算法进行求解。绕丝表面、燃料棒表面以及外套管边界均设置为绝热壁面。

因为本研究只关注流动特性，所以冷却剂定为恒物性。表1^[16]列出了LBE在350 °C时的热物性参数，该温度在CiADS的设计运行温度范围内。

图4显示了61棒束通道结构化网格无关性验证结果，采用图示a-a直线处棒束间隙的横向速度为无关性验证参数，计算工况为雷诺数等于15 000。可以看到：第3套网格(数量3 100万)与第4套网格(数量4 300万)的结果相当接近。从1 800万到4 300万四种网格两两相邻的平均横向速度偏差为0.014 0，0.014 2和0.001 0 m/s，两两网格间最大横向速度偏差为0.037 7，0.026 2和0.022 4 m/s，可见3 100万网格数量下网格与4 300万网格横向速度较前两者间已非常接近，因此选用3 100万网格进行后续计算。后续章节中Merzari验证网格以及多面体网格均经过相同步骤的无关性验证。

Merzari等^[5]使用LES方法进行了7棒束通道的流动特性研究，该结果被很多研究人员用于CFD基准题验证。在Merzari的研究中，无量纲横向流量的表达式如式(1)所示，其中z轴方向为主流方向，x轴、y轴方向以及绕丝与x轴夹角α示意图如图5显示。

式中：g_normalized(z)为轴向无量纲流量；g(z)为轴向单位高度横流量，单位为m²/s; $\overline u (z) $为轴向轴流速度，单位为m/s；l_bulk为不考虑绕丝的间隙宽度，单位为m。

定义绝对误差计算公式为

式中：g_{normalized-calulation}(z)为单位高度上截面横流量CFD计算值；g_{normalized-ANL}(z)为所参考大涡模拟计算得到的单位高度上截面横流量值，H为螺距，单位为mm。

图6显示了在直线子通道间隙处沿轴向无量纲横向流量分布对比。图中横坐标z/D为轴向位置(z)与燃料棒直径(D)的比值。可以看到：与LES模型均一化横流量与z/D积分值绝对误差最小的三项分别为：结构化网格(0.127 4)、NGR-2(0.177 2)、多面体网格(0.233 2)。多面体非结构化网格计算结果中的曲线不够平滑，这是由非结构化网格空间结构错落、不均匀所致。

图7显示了结构化网格、多面体网格的计算得到的摩擦因子与经典摩擦因子经验公式^[17-20]的对照。为方便对照，图中以结构化网格的计算结果为基准添加了±10%的误差线。结果显示：结构化网格的结果计算的摩擦因子与经验公式吻合较好，与Upgraded Cheng and Todreas(UCTD)^[17]模型的误差最小。当Re大于7 000时，二者的相对误差小于10%。同时，多面体网格得到的摩擦压降因子在Re大于7 000时逐渐偏离了经验公式预测值。这是因为相较于结构化网格，多面体网格在计算湍流时得到的涡流较弱，使得流道内轴向压降偏小。

燃料组件横截面速度流线可清晰地描述横向速度的分布情况。图8展示了61棒束通道内Z=90 mm截面上的横向速度流线分布。可以看到：两种网格的横向速度分布大体相似，但数值和位置略有偏差：(1) 结构化网格在该截面上的横向速度幅值接近0.4 m/s，而非结构化网格横向速度幅值仅为0.36 m/s，并且涡流整体大小以及分布密度上结构化网格大于非结构化多面体网格，这是造成结构化网格压力梯度计算结果大于非结构化网格的原因；(2) 结构化网格的速度幅值位于最外圈棒束的绕丝后方，而非结构化网格中位于绕丝与另一燃料棒之间的间隙上。此外，结构化网格在边通道、角通道内横向速度的涡流形态、数量都与非结构化网格不同，这是由于结构化网格划分精细且规律，可以更好地捕捉到边通道、角通道处受绕丝结构影响产生的速度涡流，单位螺距内压力梯度更大，更符合经验公式预测值。

图9显示了61棒束通道内不同子通道间的间隙编号。以中心燃料棒为第一圈，向外分别为第二到第五圈，可以将61棒束子通道间隙分为3类：I类间隙(如Ia1-Ia4、Ib1-Ib4和Ic1-Ic4)为不同圈数燃料棒间的间隙；II类间隙(如IId1-IId4、IIe3-IIe4、IIf3-IIf4和IIg4)为同一圈燃料棒之间的间隙；III类间隙(如IIIh1-IIIh5)为最外圈燃料棒与外套管之间的间隙。同时，文中参考子通道的三种类型(中心通道、边通道和角通道^[21])的划分方式进行辅助分析。

在提取不同子通道间隙的横向速度时，定义归一化横向速度U_t为

式中：l为对应该轴向位置处间隙的宽度，单位为m；u为速度，单位为m/s；$\overline u $为平均轴向速度，单位为m/s。可知，l会随着绕丝的角度不同而发生变化。

图10显示了I类间隙的无量纲横向速度U_t随轴向位置的变化关系。可以看到：同一个角度(如Ia1、Ia2、Ia3)的子通道间隙分布特性几乎重合，不同角度子通道间隙之间横向速度分布特性曲线的相位相差60°，说明横流特性在I类间隙中具有较好的周期性。这是因为中心通道的几何结构较为相似，且受外套管内壁面的影响较小。

图11显示，I类间隙Ic和II类靠近中心间隙IId1、IId2的横流特性非常相似，且越接近中心位置，相似程度越高。这是由于间隙的流动受附近燃料棒几何结构的影响较大，而受外套管内边界的影响可以忽略。结合图10，由于位于组件中心的I类(Ia1、Ia2、Ic1)、II类间隙(IId1、IId2)流动特性几乎完全一样，因此中心通道间隙特性均可用新组合公式描述；这也说明中心通道的横流特性在径向上具有较好的周期性，因而可以使用少棒束通道描述多棒束甚至无限棒束通道内的流动特性。同时看到，II类间隙最外层的IId4、IIe4间隙上的横向速度曲线较为紊乱，偏离余弦函数分布，绕丝附近横向速度的峰值大于中心通道间隙的峰值。这是由于IId4、IIe4间隙为边通道与中心通道的交界面，在边通道内形成的涡流影响范围、涡流强度都大于中心通道内的，因而绕丝附近横向速度变化更加剧烈。

图12(a)显示了III类间隙(IIIh1-IIIh5)流动特性。可以看到：III类间隙的横向速度幅值出现了明显的速度反逆转情况，这是由于当绕丝经过IIIh1间隙时，来流侧在组件边界附近形成一股比较强烈的涡流；此外，边角通道间隙横流特性曲线只有一道横向速度峰，相邻的几道边通道横流分布也较为紊乱、相似程度不高，这是由于边角通道仅受一根绕丝影响，并且由于边角通道靠着组件内边界，不同间隙附近的几何结构并非完全对称。在五组边通道间隙横流特性曲线中，由于IIIh1、IIIh5间隙临近组件边界，总体横流特性曲线也有别于相互较为相似的IIIh2-IIIh4边通道中间间隙。同时，图12(a)显示，边通道间隙横流最大值、最小值与平均值相差在20%以内，并且大部分范围内于平均值变化幅度在10%以内，基于此，本研究推荐拟合等值公式来描述边通道的横流特性。而在角通道，由于特征曲线紊乱，不容易使用方程拟合其曲线。

图12(b)显示了III类间隙在不同轴向高度上速度矢量的分布情况。同一轴向高度上，不同棒束附近的流动互有差异，在角通道附近，涡流的形态更为复杂，因而难以对边通道和角通道间隙的横流特性分布进行定量拟合。

有研究人员对带绕丝燃料棒束通道的横流特性进行了定量描述。其中，对于两中心通道间间隙横流，Naveen等^[22]提出公式：

而Wang等^[23]提出的公式为

其中：

式中：θ为绕丝相对竖直方向的倾斜角度，单位为°；α为绕丝角度，单位为°；P为燃料棒间距，单位为m；D_r为燃料棒直径，单位为m；D_w为绕丝直径，单位为m。

两位研究者公式中的u_t^*(α)定义相同，但是，Naveen等^[22]的公式没有考虑燃料棒间距(P)与燃料棒直径比(P/D)的影响，并且所采用的非结构化网格模型较为粗糙，难以得到精确结果；Wang等^[23]公式没有考虑绕丝与燃料棒间隙处的横向速度峰值。本研究对两位作者的公式进行重新整理，得到新组合公式如下：

式中：曲线的相位与绕丝位置有关，式中绕丝角度为120°；常数C取−0.208 3。

以上不同公式的计算结果与CFD模拟结果对比如图13所示。图中使用与LES基准题误差最小的结构化网格的计算结果作为以上不同公式的对照基准。可以看到：新组合公式的预测值与CFD计算值最接近，同时，能较准确地捕捉到在绕丝处的尖峰。

同时，Naveen等^[22]公式也拟合得到了边通道间的横流特性的公式：

式中：C₀=−0.152 8，C₁=4.787×10⁻⁴，C₂=−4.611×10⁻⁶，C₃=8.901×10⁻⁹。

Naveen等^[22]公式拟合自较为粗糙的非结构网格计算结果，流场信息不够细腻。本研究根据更加准确的结构化网格数据重新拟合，得到新拟合公式如下：

图14显示了Naveen公式、新拟合公式与边通道横流特性曲线与之间的对比。结果显示：相比于Naveen公式，新拟合公式更接近CFD计算结果，归一化横向速度在绕丝范围外中心值上下浮动约为20%，同时拟合式与计算值横向速度积分偏差不到1%，因此新拟合公式带来边通道的偏差对于中心间隙的流动影响很小，可认为对流动结果的影响可接受。同时可以看到，根据非结构化网格拟合的Naveen横流曲线与结构化网格模拟在边通道横流特性时相差较大。

本研究基于LES模拟基准题评估了结构化网格与非结构化网格在带绕丝燃料组件通道内流动特性的适用性，分析了网格类型对流动特性的影响，并研究了CiADS的61棒束通道内的横流特性，得到以下主要结论：

(1) 相比于非结构化的多面体网格，结构化网格的结果与LES基准题和经验公式预测值的结果误差更小。在摩擦因子计算中，结构化网格计算结果与UCTD^[17]经验公式预测值非常接近，而非结构化网格计算值偏差较大。这是由于结构化网格能够划分精细且规律，能够展示更多的横向速度涡流。

(2) 相比于非结构化网格，结构化网格CFD计算得到的横向速度幅值更大、产生的涡流更多、涡流细节更丰富。间隙横流特性主要受附近几何结构影响，表现出高度的相似性，横流分布在组件边界影响较大的间隙并逐渐偏离了正弦分布；在边通道中，无绕丝区域横向速度最大最小值与平均值相差在20%以内；而角通道间隙的横向特性呈现出强烈的无序性，角通道间隙与边通道中间间隙分布差异很大。

(3) 基于结构化网格CFD计算结果，改进了带绕丝燃料棒束中心通道与边通道横流间隙经验公式。相比Naveen公式和Wang公式，改进后的经验公式与CFD计算结果间的误差更小，有助于更准确地预测带绕丝燃料棒束横流特性。

本研究提供了子通道程序优化的新思路，同时可为CiADS燃料组件热工水力设计与分析提供参考。