液态铅铋合金在绕丝燃料棒组件中的水力特性数值模拟

姜文殊，梁峻铭，张会勇，贠相羽，卢 涛

(1.北京化工大学机电工程学院，北京 100029；2.中广核研究院有限公司暨广东省核电安全企业重点实验室，广东 深圳 518000)

目前核反应堆中应用较为广泛的金属冷却剂主要为钠的碱性液态金属、钠钾合金、锂、熔融铅和液态铅铋合金等[1].其中，液态铅铋合金由于具有熔点低、沸点高、导热性能好等一系列优点，成为了核反应堆冷却剂材料[2，3].铅基核反应堆燃料棒组件一般呈三角形排列于六边形通道中，燃料棒上含有螺旋绕丝用于燃料棒的固定.螺旋绕丝结构除为燃料棒提供结构上的支持外，还能促进流体间的横向流动以提高对流换热系数，但这种沿着加热管的螺旋形流动会增加摩擦压力损失[4].针对此类摩擦压力损失造成的压降，已有一些研究学者对此进行了相应的实验研究.

Engel[5]等对水、钠和空气在六边形燃料棒组件中的流动进行了测试，分别总结了含绕丝燃料棒束在层流、湍流以及过渡流情况下的压降，并提出了相应的摩擦系数关联式，但Engel等提出的压降关联式仅与雷诺数有关，未涉及到绕丝直径、绕丝螺距以及燃料棒直径等相关参数，因此预测偏差可能较大.

Novendstern[6]基于液态钠在217绕丝燃料棒组件中的流动实验，提出了预测绕丝燃料棒组件中压降的半经验模型，但此模型是基于绕丝均匀分布在子通道中假设，对原模型进行了一定程度的简化，可能会造成一定偏差.

Cheng[7]等同时考虑了流型和几何形状的影响，基于子通道和混合参数模型，并结合文献实验数据以及37燃料棒束流动实验数据，给出了燃料棒组件在不同流动状态下的压降关联式.

Chun[8]等以水为流动介质，在19燃料棒组件中进行了实验，并对燃料棒组件内的子通道压降进行了分析，并结合实验数据与已有的一些数据来共同评估了现有的压降关联式.结果表明，Cheng关联式与实验数据的一致性较好.

近年来，随着计算机计算及存储等能力的提升，越来越多学者开始利用CFD对绕丝燃料棒束的三维流动进行模拟.如：侯荣斌[9]等通过数值模拟主要研究了组件的压降、横向流、子通道内的流量分配、各子通道内的温度分布及燃料包壳轴向的温度分布.蔡容[10]等通过数值模拟获得了三角形排列棒束内的湍流交混系数，并拟合了新的湍流交混关系式；Natesan[11]等利用CFD数值分析计算了液态钠在19绕丝燃料组件中流动的压力损失和对流换热系数；Gajapathy[4，12]等已经通过数值模拟分别研究了P/D=1.191、1.255和1.31时，钠在217绕丝燃料棒组件中流动时的对流传热.

综上所述，目前已有的研究大多是基于水、钠等流体对燃料棒组件中流体流动的相关研究.液态铅铋合金作为目前最具发展前景的反应堆冷却剂，研究其在绕丝燃料棒组件中的水力特性十分必要.但由于液态铅铋合金的高密度、强腐蚀及实验温度要求较高等因素，开展铅铋合金在燃料组件内的热工水力实验难度较大.综合实验技术、设备及经费等原因，目前国内外较大规模地开展液态铅铋合金在燃料组件内的热工水力实验的科研机构有德国的KIT[13]和意大利的ENEA实验[14，15].为分析上述关联式对液态铅铋合金在绕丝燃料组件内流动的适用性，本文采用数值模拟方法对液态铅铋合金在19绕丝燃料棒组件中的速度分布及压降进行研究，并将模拟所得压降与Engel、Cheng以及Novendstern等压降关联式的计算结果进行对比，以验证本研究结果的准确性和可靠性，为进一步探究液态铅铋合金在绕丝燃料棒组件内的流动特性及规律积累理论基础.

1 计算模型及边界条件

1.1 几何模型

如图1所示，本文研究对象为含绕丝的19燃料棒组件，且参考了王俊[16]等给出的燃料棒直径，相关参数如表1所示.为便于分析液态铅铋合金在不同流动区域的流动特性，如图2所示，将流动区域按照不同的形状特征分为角子通道、中心子通道和边子通道.狭缝a、b、c分别表示中心子通道狭缝、边子通道狭缝及角子通道狭缝.

表1 模型参数

1.2 网格划分及网格敏感性分

由于几何模型复杂，本研究采用非结构网格划分.另外，由于绕丝与燃料棒为线接触，不利于网格划分.因此对绕丝与燃料棒进行相交简化处理，即将绕丝中心向燃料棒中心缩进0.35 mm.

为保证数值模拟计算结果的精确性和网格无关性，分别选用网格数量约为1 100万(网格1)、1 500万 (网格2)和2 100万 (网格3)的三种模型进行网格无关性验证.以狭缝a的速度分布作为网格无关性验证的基准，并对沿狭缝a的方向进行归一化处理，归一化处理后的位置用s表示，狭缝a位置如图2所示.结果表明(见图4)：三种网格所计算的沿狭缝a处的速度分布基本一致，综合考虑计算精度与计算资源，在计算精度基本相同的情况下选择网格数量较少的网格1开展后续模拟计算，网格1示意图如图3所示.

1.3 边界条件

Natesan等的研究表明，k-ε、k-ω以及Reynolds三种湍流模型所预测的摩擦因子几乎相同，对于流量为12 m3/h的湍流工况，入口雷诺数约为61 500，本文选择了standardk-ω模型进行数值计算；而流量为0.1 m3/h下的工况，入口雷诺数约为400，属于层流，因此选择laminar模型进行计算.王婧捷[17]等对不同湍流普朗特数的适用性进行了研究，结果表明，Cheng湍流普朗特数的模拟结果与实验值最为接近.因此本文选择了Cheng湍流普朗特数模型进行计算，Cheng湍流普朗特数的表达式如下.

(1)

本文不考虑燃料棒的发热，加热棒功率为0.进出口温度均设置为恒温350 ℃，并对该温度下的液态铅铋合金的物性参数进行设置，350 ℃时铅铋合金的物性参数如表2所示.入口采用速度入口，出口采用压力出口，其余壁面均设置为无滑移绝热壁面.

表2 液态铅铋合金物性参数表

2 计算结果与分析

2.1 流体流动速度

已有研究为减少网格数量，节省计算量，大多是通过周期性边界条件来对不同螺距下的流动状态进行分析，这容易受到诸如回流等出口效应的影响，导致部分真实流动特性丧失.本文采用多螺距模型对液态铅铋合金在19绕丝燃料棒组件中的流动特性进行分析，并分别提取了层流以及湍流状态下不同螺距中部和尾部截面的速度云图，如图5(层流)和图6(湍流)所示.

从图5和图6中可以看出，由于绕丝螺旋缠绕于燃料棒外表面，使得不同轴向截面位置绕丝的方位角不同，导致流动速度也会跟随着绕丝的转动方向发生相应的变化，对流体的流动起到了横向搅混作用.当六边形流道内流体流动趋于稳定后，速度云图呈现出一种周期性旋转的趋势，且绕丝方位相同的截面速度云图基本一致，这也与流道结构的周期性特点有关.比较不同子通道间的流动云图可以看出，边通道的速度明显大于中心子通道和角子通道的流速，这可能是由于边子通道的水力直径较大[18]，摩擦阻力较小，使得流动速度较大.另外，绕丝转向侧(位置1)能观察到明显的流动停滞现象，因此绕丝转向侧较另一侧(位置2)的速度不均匀性更大.

为定量分析不同子通道沿轴向方向的流体流动速度波动，分别提取层流和湍流状态下的中心子通道、角子通道和边子通道的平均速度，如图7和图8所示.从图7和图8中可以看出，第二段螺距以后，中心子通道和边子通道流速基本呈周期性变化，角子通道由于其二次流形态较为复杂，各周期速度分布略有差异，但基本是在一定范围内波动.此外，从图7和图8中可以看出，无论是层流流动还是湍流流动，边子通道沿轴向的速度均大于角子通道和中心子通道.但层流流动时角子通道沿轴向流动速度最小，而湍流流动时大概从第二段螺距开始，中心子通道的沿轴向的流动速度最小.

2.2 湍动能分析

上述湍流工况下不同螺距横截面上的湍动能分布云图如图9所示.可以看出，湍动能的分布云图与速度云图一样，均随螺距的变化呈现出一种周期性旋转的趋势，且绕丝方位相同的截面湍动能的分布基本一致.对比速度云图与湍动能分布云图后可以发现，靠近壁面处流动速度虽小，但湍动能反而更大.这是由于湍动能除了与速度有关，还与雷诺数有关，雷诺数越小，湍动能越大.而雷诺数表征的是流体惯性力与粘性力之间的比值，若把流体截面等效为无数个并联小流道，靠近壁面处流体几乎处于流动停滞状态，其雷诺数量级非常小，最后反而使得其湍动能更大.为定量分析不同子通道内湍动能沿轴向的变化情况，分别提取了不同子通道内沿轴向位置的平均湍动能，如图10所示，不同子通道内的平均湍动能与平均速度的升降趋势保持一致，且角子通道的平均湍动能要大于边子通道和中心子通道，可见，流道形状也对湍动能有着极大的影响.

图9 不同螺距下的湍动能分布图

图10 子通道平均湍动能

2.2 压降分析

本文将模拟所得到的压降结果与Engel关系式、Cheng关系式以及Novenstern关系式的计算结果进行对比，以验证本研究结果的准确性和可靠性，同时评估不同压降关系式对液态铅铋合金绕丝燃料组件的适用性.

对于不可压缩流体，压降系数的计算公式为

(2)

公式中：L为压降测量截面的距离；De为水力直径，ν为流体流速；f为摩擦因子.

对于摩擦因子的计算，Engel[5]关系式为

层流时

(3)

湍流时

(4)

过渡流动时

(5)

公式中：Ψ为间断因子，是一个连接湍流区和层流区的参数，Ψ=(Re-400)/4 600.

Cheng[7]关系式为

层流(RebL≥Reb)时

(6)

湍流(RebT≤Reb)时

(7)

过渡流(RebL

(8)

在上述Cheng关系式中

(9)

(10)

(11)

(12)

(13)

以上诸公式中：P为燃料棒间距；D为燃料棒直径；H为绕丝螺距.

Novenstern关系式

fi=fsM，

(14)

公式中：fi为子通道i的摩擦因子；fs为光滑管道的摩擦因子；M为乘法因子，是基于光管摩擦因子的修正系数，它主要受P/D和H/D的影响，其计算公式为

(15)

若光滑管道中流动的摩擦因子fs通过式(16)进行计算

(16)

则压降计算公式为

(17)

公式中：VT为平均流体流速；X为流动分布因子，它可以通过中心子通道、边通道以及角通道计算，其计算公式为

(18)

公式中：Ni和Ai分别是子通道i的数量和面积；Dei为子通道i的等效直径.

若采用Blasius关联式进行近似，光滑管道中流动的摩擦因子为

(19)

则压降公式为

(20)

本文分别提取了层流和湍流状态下角子通道，边子通道以及中心子通道的压降的变化，如图11所示.从图11中可以看出，无论是层流还是湍流状态下，角子通道，边子通道以及中心子通道的压降的变化趋势基本一致.不同的是，层流条件下压降变化较小，且基本呈线性变化；这可能是由于层流速度较小，层流状态下相同距离截面间受绕丝位置的影响不明显，因此压降呈线性变化.而湍流状态下的流体压降变化很大，且随螺距呈现一种周期性变化见图11(b).这是由于湍流速度大，绕丝不同位置对流体流动的速度影响不同，导致相同截面间的压降存在显著差异，因此湍流状态下压降呈现为周期性变化.

图11 子通道压降对比

为验证本研究压降的准确性和可靠性，选取中心子通道平均压降分别与Engel关联式、Cheng关联式以及Novendstern关联式进行对比，如图12和图13所示.

图12 层流子通道压降关联式对比

图13 湍流子通道压降关联式对比

从图12和图13中可以看出，层流状态下，本文压降模拟值与Novendstern关联式预测的压降结果更为接近；湍流状态下，压降的模拟值与Cheng关联式和Novendstern所计算的结果符合度比较高.湍流状态下，因为Engel关联式属于纯经验公式，未考虑绕丝螺距、燃料棒直径以及燃料棒间距等因素的影响，其准确性较其他关联式较低.Novendstern弥补了Engel关联式和Cheng关联式的不足，不仅考虑了绕丝螺距、燃料棒直径以及燃料棒间距等因素的影响，还考虑了绕丝导致的流体横向分布的影响，因此无论层流还是湍流状态下，Novendstern关联式均能较为准确的预测压降.

3 结论

本文利用数值模拟方法对液态铅铋合金在燃料棒组件内流动进行了数值模拟，并分别对液态铅铋合金在绕丝燃料棒组件中的速度分布及轴向压降进行了研究，并与已有的压降关联式进行了对比，可以得到以下结论：

(1)从第二段螺距开始，速度分布基本呈现一种周期性变化.子通道速度分析表明，无论是层流流动还是湍流流动，边子通道的平均速度均大于角子通道和中心子通道.但层流流动时角子通道平均速度最小，而湍流流动时中心子通道平均速度最小.

(2)无论是层流还是湍流状态下，角子通道，边子通道以及中心子通道的压降的变化趋势基本一致.不同的是，层流压降变化很小，且基本呈线性变化；而湍流状态下的流体压降变化很大，且随螺距呈现一种周期性变化.

(3)压降数值模拟结果与已有关联式计算结果对比表明，本文的研究结果具有较高的准确性.层流状态下，模拟值与Novendstern关联式预测的压降结果更为接近；湍流状态下模拟值与Cheng关联式和Novendstern关联式的计算结果的符合度较高.