铅铋冷却快堆含绕丝燃料组件子通道程序开发与验证

刘佳泰，彭天骥,3，苏兴康，顾 龙,4,*

(1.中国科学院 近代物理研究所，甘肃 兰州 730000；2.中国科学院大学，北京 100049；3.先进能源科学与技术广东省实验室，广东 惠州 516003；4.兰州大学，甘肃 兰州 730000)

中国加速器驱动嬗变研究装置(CiADS)由直线加速器、散裂靶和次临界反应堆组成[1]，其中次临界反应堆选用液态铅铋(LBE)冷却快堆[2]。铅铋冷却快堆具有较高的安全性和可持续性等优点，然而受限于LBE非线性的湍流普朗特数以及LBE具有腐蚀性等困难，通过实验获得堆芯的精细热工水力参数所需要满足的条件十分严格。CFD的几何刻画精细，但对整个堆芯或燃料组件进行网格划分和数值计算需要大量的计算资源，计算时间较长。子通道分析方法通过将燃料组件划分为内通道、边通道和角通道，将经验关系式代入单个子通道的守恒方程中使方程组封闭，迭代求解获得一定精度的热工水力参数，可大幅减少精细计算时间。

对于液态金属，国内外开发了许多适用于含绕丝燃料组件的子通道程序。早期多是基于钠冷快堆，如ASFRE-Ⅲ[3]基于分布阻力模型对钠冷快堆进行计算；MATRA-LMR[4]是基于压水堆程序MATRA开发的钠冷快堆子通道程序；COBRA-WC[5]和COBRA-LM[6]均是基于美国西北太平洋实验室开发的COBRA程序改进成的钠冷快堆程序。近年来针对铅铋冷却快堆的研究日益增多，如SACOS-PB[7]可计算铅冷快堆稳态条件下的温度分布；ANTEO+[8]是一种适用于液态金属强制对流的通用程序。目前已有的液态金属子通道程序由于使用的关系式不同，具有不同的适用条件和限制，是否适用于CiADS的子通道分析需进一步验证。

本文开发适用于铅铋冷却快堆的子通道程序，对液态铅铋的摩擦阻力模型、湍流交混模型和对流换热模型进行适用性分析，并与含绕丝燃料组件的LBE大涡模拟(LES)计算结果和传热实验结果进行对比验证。

1 守恒方程及本构关系

本研究以铅铋冷却快堆含绕丝燃料组件为研究对象，采用子通道分析程序对稳态强迫对流工况下冷却剂及包壳表面的温度分布进行计算。程序基于Fortran平台进行开发，由于LBE在正常工况下不会产生沸腾，因此采用单相流状态下的守恒方程。

1.1 质量守恒方程

(1)

式中：A为轴向流动面积；ρ为冷却剂密度；t为时间；m为轴向质量流量；z为轴向控制体高度；eik为符号函数，i为子通道的编号，k为子通道i周边的间隙；w为单位长度横向质量流量。方程各项分别为质量随时间的变化、轴向质量流量随空间的变化和间隙处的横向质量流量。

1.2 轴向动量守恒方程

(2)

式中：U为轴向流速；p为压力；g为重力加速度；θ为冷却剂通道与垂直方向的夹角；f为摩擦阻力系数；D为燃料棒直径；K为形阻系数；fT为湍流动量交换系数；w′为单位长度湍流交混流量。方程左边3项分别为轴向动量随时间的变化、轴向动量通量随空间的变化和横流导致的横向动量变化。方程右边分别为压力项、重力项、摩擦力项和湍流交混项。方程右边各项共同形成轴向力。

1.3 横向动量守恒方程

(3)

式中：s为间隙宽度；l为子通道i与子通道k之间的质心距离；KG为横向阻力系数；Δp为相邻子通道间的横向压差。方程左边分别为横向流量随时间的变化和横向流量通量随空间的变化，方程右边分别为横向压差项和间隙阻力项。由于绕丝结构产生横流，使燃料组件的同一横截面内具有明显的压力差。

1.4 能量守恒方程

(4)

式中：h为流体焓值；PW为燃料棒r面向子通道i的加热周长；φir为燃料棒r面向子通道i的份额；q为燃料棒的线功率密度；CQ为冷却剂产热份额；φin为燃料棒与通道之间的接触份额；q′为燃料棒传入冷却剂的线功率密度；Δh为相邻通道间焓差；Ck为横向换热系数；ΔT为相邻通道间的温度差。方程左边分别为焓随时间的变化、轴向焓通量随空间的变化和通道内所有间隙处的横向焓通量。方程右边分别为燃料棒传入流体的热量、直接在冷却剂中产生的热量、湍流交混的热量和相邻通道间的热量交换。

1.5 计算流程

子通道程序在轴向节点的入口沿流向逐步计算。首先通过能量方程更新焓值，之后对质量方程和动量方程联立的矩阵采用高斯-赛德尔迭代计算得到横向流量。最后将横向流量带入到质量方程计算得到轴向流量，带入到横向动量方程计算得到压力梯度。程序求解流程如图1所示，主程序在读取边界条件和几何条件后，依次调用物性模块、能量方程迭代模块和动量方程迭代模块完成循环，直到轴向节点全部计算收敛后输出计算结果。

图1 子通道程序中冷却剂求解流程

1.6 物性模型

对程序模型进行不确定性分析是评估和验证模型准确性的关键环节。冷却剂的物性作为模型关键输入参数，在不确定性的传播过程中起到关键作用。冷却剂相关物性的不确定度列于表1[9]，包括密度、比定压热容、动力黏度及热导率。程序中还需已知冷却剂的焓值，通过焓的定义计算得到。

表1 LBE的推荐物性关系式

(5)

式中：h0为熔点温度处的焓值；T为冷却剂温度；TM为熔点；cp为比定压热容。

1.7 摩擦阻力模型

程序采用达西公式计算压降：

(6)

式中：L为通道长度；De为通道水力直径；v为冷却剂流速。

求解压降的关键是得到合适的摩擦阻力系数f，而摩擦阻力系数一般是雷诺数Re和几何关系的函数。Fan等[10]以CiADS棒束为原型，开展了阻力系数测量实验，并将实验结果与摩擦关系式进行对比。如图2所示，在3 750≤Re≤16 250范围内，Rehme关系式[11]能在±10%的范围对全部54个实验点准确预测，均方根误差为7.8%。因此本程序选用Rehme关系式计算摩擦阻力系数：

图2 Rehme关系式与CiADS实验的摩擦阻力系数对比

(7)

式中：Pwb为棒束湿周；Pwt为总湿周；F为几何因子，用于考虑绕丝的几何结构。

(8)

式中：P为相邻燃料棒质心的距离；Dw为绕丝直径；H为绕丝螺距。关系式的适用范围为：103≤Re≤3×105，1.125≤P/D≤1.417，5≤H/D≤50。

1.8 湍流交混模型

绕丝的存在会增加相邻子通道间的横向交混，具有展平横向温度分布的作用。假设湍流交混不引起质量交换，只引起动量和能量上的交换。考虑湍流交混时的横向流量，采用以下公式：

(9)

内部通道为：

(10)

(11)

(12)

壁面处通道为：

C1L,Lam=0.413(H/D)-0.5(Ar2/A′2)0.5tanθ

(13)

C1L,Tur=0.73(H/D)-0.5(Ar2/A′2)0.5tanθ

(14)

(15)

其中：

ψb=lg(Reb/RebL)/lg(RebT/RebL)

(16)

(17)

(18)

图3 湍流交混关系式计算的交混系数与实验结果对比

1.9 对流换热模型

对于包壳外表面和冷却剂之间的对流换热，需要通过努塞尔数(Nu)来确定对流换热系数。对于LBE，Nu通常被认为是棒径比和贝克莱数(Pe)的函数。Pacio等[13]在LBE传热实验中选取了3个加热段测量位置ML1(z/H=1/6)、ML2(z/H=11/6)和ML3(z/H=15/6)，并与液态金属的对流换热关系式进行对比，最终推荐Kazimi等[14]关系式用于计算液态铅铋的努塞尔数，如式(19)所示。图4给出经验关系式与实验测得努塞尔数之间的对比，可看出关系式低估了ML1处的Nu，这是因为ML1在热发展区，但在位于充分发展区的ML2和ML3时对Nu的预测比较准确，均方根误差控制在7.1%。因此程序采用Kazimi等关系式计算努塞尔数，适用范围为：10≤Pe≤5 000，1.1≤P/D≤1.40。

图4 对流换热关系式计算的努塞尔数与实验结果对比

(19)

2 程序验证

2.1 子通道程序流动验证

采用CFD可获得精细的流场数据，为验证CFD结果可靠，采用大涡模拟计算结果作为基准验证[15]，该结果为美国阿贡国家实验室通过开源程序Nek5000计算得到，具有高保真度。本研究在此模型基础上进一步划分子通道供程序验证使用。

棒束结构及子通道编号如图5所示，绕丝沿流向逆时针方向旋转。表2列出棒束的具体几何参数。几何建模时，由于绕丝与燃料棒之间为线接触，不利于网格划分，因此在每个绕丝与燃料棒接触处做0.25 mm的倒圆角处理，以提高网格质量。

图5 7棒束子通道编号及划分

表2 7棒束几何参数

采用商业计算流体软件STARCCM+进行多面体网格划分。图6示出整体网格及间隙处的局部网格细节。在近壁面处设置4层边界层网格，选取一个螺距高度进行计算，进出口设置周期性边界，进口雷诺数为9 457，生成网格后，使用SSTk-ω湍流模型进行计算。

图6 7棒束网格划分上视图

为比较结果，将间隙横流速度u在整个流向上进行积分，表示为坐标s的函数：

(20)

分别选取0.8 mm(网格1)、0.6 mm(网格2)、0.3 mm(网格3)和0.2 mm(网格4) 4套不同尺寸的网格进行无关性验证。图7为网格计算结果与LES结果(ANL-LES)在A-A、B-B、C-C间隙处的横向速度对比。网格3与网格4的结果相近且与LES结果吻合，证明基于大涡模拟流场数据的CFD模拟结果较好，可基于网格3划分子通道进一步验证子通道程序。

a——间隙A-A；b——间隙B-B；c——间隙C-C

为验证程序流动计算的准确性，将子通道程序的计算结果与CFD结果进行对比。图8为出口处不同子通道的质量流量分布，边通道的出口流量最大，内通道其次，角通道出口流量最小，程序能准确计算出口质量流量。图9为子通道1的冷却剂轴向速度沿流向在1个螺距长度的分布，可看出子通道程序计算的轴向速度与CFD计算结果具有相同的趋势。

图8 出口处子通道质量流量分布

图9 子通道1轴向速度分布

2.2 子通道程序传热验证

子通道程序要求能准确计算冷却剂温度和燃料棒包壳外表面温度。为验证程序传热计算的准确性，选取Pacio等[13]开展的19棒束含绕丝燃料组件的棒束传热实验进行验证。燃料组件的具体几何参数列于表3。具体边界条件包括进口温度Tin=473 K、进口质量流量M=19.18 kg/s及总加热功率Q=197.04 kW，加热功率均匀分布在每根棒的加热段上。在程序中模拟824 mm的流动发展段和870 mm的加热段，棒束子通道划分与编号方式如图10所示。

图10 19棒束的子通道编号和划分

包壳和冷却剂温度计算结果与实验结果对比如图11～13所示，分别对应沿加热段轴向高度的3个不同测量位置ML1、ML2和ML3。其中ML1处于热发展段，ML2和ML3处于热平衡段。除实验数据外，另外选取子通道程序SACOS-PB的计算结果进行对比[16]。由图12、13可看出，内通道比边通道和角通道更热，ML3比ML2的温度曲线变化更明显。同一截面上，边通道和角通道的温度均低于流体平均温度。程序计算结果与实验数据在子通道处的最大相对误差为4.08%，包壳外表面处的最大相对误差为4.36%，且与其他子通道程序的计算精度类似，总体计算相对误差低于5%，验证结果较好。

图11 ML1处包壳外表面和冷却剂温度分布

图12 ML2处包壳外表面和冷却剂温度分布

3 结论

本文开发了适用于液态铅铋冷却含绕丝燃料组件的子通道分析程序，并对不同流动和传热关系式进行了评价。基于7棒束大涡模拟流场数据验证了程序流动模块的有效性，基于19棒束传热实验数据验证了程序传热模块的有效性。结果表明开发的子通道程序能较好地计算液态铅铋冷却绕丝定位燃料组件的流动传热特性，可代替或辅助CFD计算，节约了计算时间。

图13 ML3处包壳外表面和冷却剂温度分布