基于四方程模型的LBE环管湍流换热数值研究

苏兴康，顾 龙,3,*

(1.中国科学院 近代物理研究所，甘肃 兰州 730000；2.中国科学院大学 核科学与技术学院，北京 100049；3.兰州大学 核科学与技术学院，甘肃 兰州 730000)

由于铅及铅铋共晶合金(LBE)热工水力性能良好、低压下沸点高、化学惰性强，铅冷快堆(LFR)被选为第4代核能系统的主力堆型之一[1]。中国加速器驱动嬗变研究装置(CiADS)将LBE冷却快堆作为次临界堆堆型[2-3]。三角形排列棒束是LFR燃料组件常见的形式，研究LBE冷却棒束的热工水力特性对LFR安全性和经济性具有重要意义。在不考虑燃料棒相互作用的情况下，三角形棒束可简化为用六边形传热单元等效的环管[4-6]。因此，环管内的湍流换热研究有助于分析三角形棒束热工水力。LBE的不透明、腐蚀性特点增加了LBE环境中传热实验的难度[7-8]，需借助计算流体力学(CFD)数值模拟LBE在环形管道内的湍流换热传热过程。然而，LBE的分子普朗特数(Pr≈0.01～0.03)较低，导致雷诺时均方法(RANS)中的雷诺比拟假设会影响LBE的传热计算[9]。即用于模拟传统流体(Pr≈1)的恒定湍流普朗特数(Prt)，可能会影响LBE传热的评估[10]。在过去的40年里，为提高LBE的计算精度，国内外研究了用于封闭能量方程的两方程kθ-εθ模型。用于封闭动量方程的两方程k-ε湍流模型和用于封闭能量方程的两方程kθ-εθ传热模型统称为四方程k-ε-kθ-εθ模型。

文献[11-21]对四方程模型的系数函数和近壁热湍流行为进行了大量的研究，发展了多种适用于不同Pr流体的四方程模型。近年来，Manservisi等[22-24]在前人研究的基础上，发展了新的四方程LBE湍流换热模型，并利用直接数值模拟(DNS)数据和LBE传热实验数据对模型进行了验证，在平板、管道、三角形和四边形棒束通道内的计算结果较好。文献[25-30]在Manservisi四方程模型的基础上，提出了比耗散率k-ω-kθ-ωθ和对数比耗散率k-Ω-kθ-Ωθ形式的四方程模型，并进一步计算和验证了LBE在平面、管道、后台阶和六边形棒束中的湍流换热，预测精度良好。

尽管人们对低Pr流体的可靠计算工具愈加感兴趣，但商业代码仍较缺乏，基于四方程模型的LBE环管内的热工水力研究较少。因此，本文基于开源CFD软件OpenFOAM，开发四方程求解器4eqnFoam。基于4eqnFoam，计算LBE在平板和圆管内的湍流换热，并与DNS数据和实验关系式进行对比验证。最后，基于该求解器，计算LBE在环管内的流动传热，获得不同加热条件下LBE的关键参量，为计算铅铋流体湍流换热提供参考。

1 湍流换热模型

环管几何模型示于图1。内、外壁面分别受恒热流qi、qo加热。边界条件划分为入口Γi、出口Γo、外壁面Γw,o和内壁面Γw,i。Ri和Ro分别为内外半径。

图1 同心环管示意图Fig.1 Sketch of concentric annulus

对于不可压缩、常物性、无重力的充分发展雷诺时均方程组，有：

(1)

(2)

(3)

(4)

其中，ux为沿流动方向的速度分量。

(5)

其中，k和νt分别为湍动能和湍流运动黏度。采用引入了Kolmogorov速度尺度uε=(vε)1/4的Abe两方程k-ε湍流模型计算νt[17]。ε为湍流耗散率。

(6)

其中，αt为湍流热扩散系数。选用Manservisi两方程kθ-εθ传热模型计算αt[22-24]。该模型引入了湍流温度动能kθ和湍流温度耗散率εθ，能较好地计算LBE湍流传热。

2 数值求解器

OpenFOAM是基于C++编写的开源CFD类库[31-32]。基于上述湍流传热模型和OpenFOAM，开发了4eqnFoam求解器。如图2所示，4eqnFoam主要由动量方程UEqn.H、能量方程TEqn.H、SIMPLE压力泊松方程pEqn.H、k-ε湍流模型VEqn.H、kθ-εθ换热模型HEqn.H等子程序组成。前处理时，采用GAMBIT软件划分网格得到mesh文件，通过fluentMeshToFoam工具导入OpenFOAM中得到网格文件polyMesh。在前处理时，初始及边界条件数据、网格及流体物性数据、离散及求解格式数据分别存放于0、constant、system文件夹中。求解器执行计算时，前处理数据被求解器相关子程序调用。基于有限体积法离散并求解稳态控制方程。对流项采用高斯迎风格式、扩散项采用高斯线性格式。矩阵选用高斯赛德尔进行迭代求解，残差收敛条件为：Max|Qn+1/Qn-1|<10-10，其中Q代表ui、T、p、k、ε、kθ、εθ，n为迭代步数。

图2 4eqnFoam求解器在OpenFOAM中的框架Fig.2 Framework of 4eqnFoam solver in OpenFOAM

3 计算结果

3.1 数值验证(三维平板)

图3 平板无量纲速度、无量纲雷诺应力、无量纲温度和无量纲雷诺热流Fig.3 Non-dimensional velocity,non-dimensional Reynolds stress,non-dimensional temperature and non-dimensional Reynolds heat flux of plane

3.2 数值验证(圆管，Ri=0)

无量纲速度u+和无量纲雷诺应力τ12沿圆管壁面无量纲法向距离y+的分布如图4所示。由图4a可见，工况B、D与管道壁面速度分布规律吻合良好，线性区满足u+=y+，且对数区满足u+=2.5lny++5，相对误差均在4%以内。由图4b可见，工况A无量纲雷诺应力分布与工况DNS-A吻合良好，相对误差在2%以内。

图4 圆管无量纲速度和无量纲雷诺应力Fig.4 Non-dimensional velocity and non-dimensional Reynolds stress of pipe

图5 圆管无量纲温度和无量纲温度脉动Fig.5 Non-dimensional temperature and non-dimensional temperature fluctuation of pipe

圆管无量纲雷诺热流〈uj′θ′〉+的分布如图6a所示，工况A的雷诺热流分布预测结果与DNS-A分布趋势相近，整体预测结果高于DNS结果，最大相对误差约为20%。努塞尔数Nu随派克利数Pe的分布示于图6b，将努塞尔数计算结果与Kirillov推荐的圆管传热关系式(Nu=4.5+0.018Pe0.8，25

3.3 环管内流动传热

1)环管Ⅰ：外壁绝热内壁恒热流加热(qi=qw，qo/qi=0)

表1 环管Ⅰ计算工况Table 1 Calculation condition of annulus Ⅰ

图7为环管Ⅰ无量纲速度u/ub和无量纲雷诺应力τ12分别沿外壁离开到内壁的无量纲距离y/d的分布。工况A的速度场、雷诺应力场均与工况DNS-A吻合良好。最大速度出现在y/d=1.10处，零雷诺应力出现在y/d=1.14附近。

图7 环管Ⅰ无量纲速度和无量纲雷诺应力Fig.7 Non-dimensional velocity and non-dimensional Reynolds stress in annulus Ⅰ

图8 不同雷诺数下的环管Ⅰ无量纲温度和无量纲温度脉动Fig.8 Non-dimensional temperature and non-dimensional temperature fluctuation for different Re in annulus Ⅰ

图9a、b分别为不同雷诺数下环管Ⅰ的无量纲雷诺热流〈uj′θ′〉+和无量纲热扩散系数α+=αt/α(α+用于表征湍流热流量与分子热流量的比值)。由图9a可见，雷诺数越大，无量纲雷诺热流峰值越大。由图9b可见，在近壁区，α+≪1，分子导热作用占主导地位。当LBE的流动雷诺数较小时，全流场α+均小于1。只有当雷诺数增加到一定程度时，如工况F，开始出现转折点α+≥1，湍流热扩散作用在主流区开始强于LBE的分子导热作用。

环管Ⅰ中的湍流普朗特数Prt=vt/αt分布如图9c所示。主流区的Prt分布相对近壁区平滑。随着雷诺数的增加，平均湍流普朗特数〈Prt〉逐渐减小，当雷诺数增加到一定程度时，如工况D、E和F，〈Prt〉变化不明显。此外，选取了多个外壁绝热内壁加热的环形管道内液态金属传热关系式对比Nu计算结果，如文献[38-42]的传热关系式。在低Pe下4eqnFoam预测的Nu结果与仇子铖[39]和Jaeger[40]关系式符合良好，中等Pe下与朱锋杰[38]和石双凯[41]关系式符合良好，而高Pe下与张贵勤[42]传热关系式符合良好(图9d)。

图9 不同Re下环管Ⅰ的无量纲雷诺热流、无量纲热扩散系数、湍流普朗特数和努塞尔数Fig.9 Non-dimensional Reynolds heat flux,non-dimensional thermal diffusivity,turbulent Prandtl number and Nusselt number for different Re in annulus Ⅰ

2)环管Ⅱ：内外壁面热流密度相等(qi=qo=qw)

令qi=qo=qw，环管Ⅱ其他计算工况与表1一致。内外壁面热流密度相等时，θ+的最大值和零梯度值出现在最大速度值对应位置(y/d=1.10)附近(图10a)。相比于外壁绝热时，无量纲温度脉动在外壁附近多出1个峰值，整体波动变化情况更加剧烈，且两峰之间的谷值位置均偏向内壁(图10b)。无量纲雷诺热流的分布曲线和雷诺应力分布曲线均不对称(图10c)。零雷诺热流的位置出现在y/d=1.07附近，且不随雷诺数变化。随着雷诺数的增加，平均湍流普朗特数〈Prt〉逐渐减小，当雷诺数增加到一定程度时，如工况D、E和F，〈Prt〉变化不明显(图10d)。

图10 不同Re下环管Ⅱ的无量纲温度、无量纲温度脉动、无量纲雷诺热流和湍流普朗特数Fig.10 Non-dimensional temperature,non-dimensional temperature fluctuation,non-dimensional Reynolds heat flux and turbulent Prandtl number for different Re in annulus Ⅱ

3.4 四方程模型求解器性能分析及适应性评价

当前，4eqnFoam求解器是基于OpenFOAM V6版本特性编写的，编译集成在Ubuntu 18.04系统环境中。串行运行下，求解不同几何和网格算例时，求解器所消耗的计算时间和内存列于表2。测试处理器为Intel(R)Core(TM)i7-9700 CPU @ 3.00 GHz，系统运行内存为16 GB @ 2667 MHz，测试计算平台是8核心8线程的HP ProDesk 680 G4 MT。

表2 4eqnFoam求解器性能分析Table 2 Performance analysis of 4eqnFoam solver

从平板、圆管和环管的计算验证结果来看，本研究开发的4eqnFoam求解器具备了初步计算不可压缩、常物性和无重力条件下LBE湍流换热过程的能力，并能获得雷诺热流、温度脉动和湍流普朗特数等结果。由于验证过程中只考虑了Pr=0.025常物性下LBE在简单几何中的对流传热，对于更广范围内的低Pr流体、复杂几何结构(后台阶、棒束通道等)以及变物性条件下4eqnFoam的适应性和计算精度，还需进一步研究。

4 结论

为获得LBE在环管内不同加热条件下的传热现象，基于开源CFD程序OpenFOAM，开发了四方程模型求解器4eqnFoam。基于4eqnFoam求解器，将LBE在平板和圆管内的传热计算结果与DNS和实验关系式进行了对比，验证了4eqnFoam的有效性。将模拟计算扩展到LBE在内壁加热外壁绝热、内外壁面等热流加热的环管内，获得了不同雷诺数下LBE充分发展的无量纲温度、温度脉动、雷诺热流、热扩散系数和湍流普朗特数等分布。结果表明：随着雷诺数的增加，两种加热条件下LBE在环管内的无量纲温度脉动、雷诺热通量、湍流热扩散作用变大，而平均湍流普朗特数减小；靠近壁面的湍流普朗特数较湍流核心处的变化剧烈；当雷诺数增加到一定程度时，平均湍流普朗特数变化不明显。基于四方程模型求解器，可为计算LBE热工水力现象提供更多参考。