Vo  l  . 55 ,No. 6
Jun.2021
第55卷第6期
2021年6月
原子能科学技术
AtomicEnergyScienceandTechno  ogy
基于OpenFOAM 的液态金属铅铋三维流动
换热特性数值模拟研究
何少鹏王明军，章静田文喜°，苏光辉12，秋穗正'
"•西安交通大学动力工程多相流国家重点实验室，陕西西安710049#
2.西安交通大学核科学与技术学院，陕西西安710049)
摘要：铅铋快堆是第4代核能系统的主要堆型之一，但由于液态金属铅铋的热物性与传统工质如水、空
气等有很大不同，假设流动边界层与热边界层相似的雷诺比拟原理已不再适用％本文在开源程序
OpenFOAM 中开发了基于
四因子模型的自定义求解器，考虑热边界层与流动边界层的差异
性，对带绕丝棒束通道中液态金属铅铋的流动换热现象进行数值模拟，得到了速度、温度等重要热工水 力参数的三维分布，揭示了绕丝对冷却剂流动传热过程的影响规律，并将计算结果与经典实验关联式进 行对比，结果符合良好，证明了所用模型和程序的正确性％本研究可为在OpenFOAM 中添加新模型、开
发自定义求解器以及开展针对液态金属流动换热问题的计算流体动力学(CFD)模拟提供参考％ 关键词：液态金属；湍流换热模型；计算流体动力学;OpenFOAM
中图分类号:TL331 文献标志码:A  文章编号：1000-6931(2021)06-1007-08
doi ：10. 7538/yzk. 2020& youxian. 0445
Numerical  Simulation  of  Three-dimensional  Flow  Heat  Transfer  Characteristics  of  Liquid  Lead-bismuth  Using  OpenFOAM
HE  Shaopeng 12, WANG  Mingjun 12*  , ZHANG  Jing'2, TIAN  Wenxi 1'2,
SU  Guan g hui 1 2 !QIU  Suizhen g 12
收稿日期：2020-06-28#修回日期：2020-08-06基金项目：国家自然科学基金资助项目(11705139)
*通信作者：王明军
(1. State  Key  Laboratory  on  Poxver  Engineering  and  MultiPhase  Flov  ,
Xi^an  Jiaotong  University  , Xi  an  710049 , China  #
2. School  of  Nuclear  Science  and  Technology  , Xi' an  J  iaotong  University  , Xi  an  710049 , China )
Abstract : Liquid  metal  reactor  is  one  of  the  main  types  of  the  Gen-' nuclear  power
system. However  , because  the  thermal  properties  of  liquid  lead-bismuth  are  very  differ -
entfromthoseoftraditionalworkingmediasuchaswaterandair  theReynoldsanalogy  principle  whichassumesthattheflowboundarylayerissimilartothethermalbounda-
rylayer  isnolongersuitable.Inthispaper  acustomsolverbasedon kI )Ik I )-four- parametermodel  wasdevelopedin  opensourceprogram  OpenFOAM.Thedi f erence  between  the  thermal  boundary  layer  and  the  flow  boundary  layer  was  considered. The
1008原子能科学技术第55卷
flow heat transfer simulation of liquid lead-bismuth in a wire-wrapped rod bundle was carriedout.Thethree-dimensionaldistribution ofsomeimportantthermalhydraulic parameterssuchasvelocityandtemperaturewasobtained.Theinfluenceofwireonthe flow and heat transfer process of coolant was revealed,and the calculated results were compared with the classical experimental correlations.The results are in goodagree­ment,which indicates the correctness of the model and program.This study may provide references for CFD simulations of liquid metals flow and heat transfer by adding new models and developing custom solvers in OpenFOAM.
Key words：liquid metal；turbulence heat transfer model；CFD；OpenFOAM
液态金属铅铋快堆作为第4代核能系统的主要堆型之一，有自然循环能力强、中子经济性好等诸多优点。但液态金属铅铋由于特殊的热物性质(如低普朗特数、高热导率等)，其流动边界层与热边界层已不在一个量级，流动换热特性已与水、空气等常见工质有很大不同[1]，这给液态金属铅铋反应堆的研究设计带来了许多困难％故亟需针对液态金属铅铋开发更为精确的湍流换热模型以进行高精度数值模拟％
在主流商用CFD软件中，通常采取定值湍流普朗特数的SED模型对液态金属流动换热问题进行模拟求解。该模型基于雷诺比拟原理，即假设热边界层与流动边界层具有相似性，其成立的条件为P&接近于1%对于空气、水等常见工质，使用湍流普朗特数为0.85〜1的SED模型往往能得到工程上可接受的结果%但针对液态金属等低普朗特数流体开展的一些CFD模拟研究表明，定值湍流普朗特数的SED 模型已不能达到很好的预测精度％Kawamura 等2利用直接数值模拟研究了雷诺数与普朗特数对液态金属流动换热的影响，结果表明采用定值湍流普朗特数的SED模型已不能很好地与实验关联式吻合％Cheng等囚通过CFD模拟研究了单通道内液态金属的热工水力现象，发现当SED模型的湍流普朗特数由0.9增加到1.5时，预测努塞尔数将会下降约30%% Chen等4使用一种新网格标记分离方法研究了绕丝棒束组件三维流动换热特性，结果表明绕丝会引起横流以及轴向速度周期振荡，影响热工水力特性％Manservisi等5针对液态金属开发了一种更为精确的k-e-k e-e e四因子模型%该模型是一种显式的代数热通量模型，充分考虑了热边界层与流动边界层的差异性，引入脉动温度方均值及其耗散率两个新参数，建立了湍流热扩散率与脉动温度方均值间的数值关系，能实现对液态金属湍流换热的精确模拟% Manservisi等时*使用四因子模型在三角形栅元通道和四边形栅元通道内进行了数值模拟与验证，均取得了较好的预测结果％
考虑到主流商业软件代码闭源，封装严格，不利于用户对模型进行修改和开发，而开源计算流体力学软件OpenFOAM采用C++语言编程，具有编程环境自由、模型开发简便等优点，广泛应用于能源化工、海洋船舶等领域％基于开源平台进行模型添加和求解器开发对实现更高精度的液态金属流动换热CFD模拟有重要意义％故本文基于h四因子模型在开源程序OpenFOAM 中开发针对液态金属的自定义求解器，阐述数值模拟中所使用的数学物理模型，包括守恒方程与四因子模型方程，并对带绕丝棒束通道中液态金属铅铋的流动换热问题进行模拟研究%
1数学物理模型
reactor types1.1守恒方程
液态金属反应堆中液态金属冷却剂在正常工况下的流动可视为单相、不可压缩的湍流流动％因此，经过不可压缩流简化及雷诺时均处理后，动量、质量、能量守恒方程如式(1)所示％
-(UU)=
▽+•(U'U') ,•U=0
3(4+3^•(U T)=
.V•(#T)—3^•(YY)
(1)
第6期何少鹏等：基于OpenFOAM的液态金属铅铋三维流动换热特性数值模拟研究1009
式中!为速度平均值；T为温度平均值;U'为湍流脉动速度；T'为湍流脉动温度；C”为比定压热容#为热导率％
除要求解的U、P、T3个未知量外，经雷诺时均处理后，方程组(1)还引入了另外2个额外未知量，即湍流雷诺应力张量！可和湍流热流失量！T"，统称为湍流脉动附加项间。为使方程组封闭可解，需补充相应的模型和方程来求解2个新增的未知量％对于湍流雷诺应力张量，一般采用湍流模型来求解，在本文中，SED 模型应用标准k-模型，四因子模型则应用低雷诺数修正的k-模型％对于湍流热流矢量，将其表示为湍流热流矢量输运方程的代数近似解的形式，即：
YY=—&P T(2)经由式(2)，未知量湍流热流矢量！T被替换为湍流热扩散系数&%对湍流热扩散系数的计算，是四因子模型和SED模型的主要不同之处％在基于雷诺比拟假设的SED模型中，湍流热扩散系数按式(3)计算％
Pr t=v t/a t(3)湍流普朗特数给定后(本文分别取针对液态金属推荐的湍流普朗特数15和2.0进行模拟计算)，借由湍流模型求得的局部湍流黏度便可进一步求解出湍流热扩散系数至此，方程组(1)得以封闭，方程组理论可解％
1.2四因子模型
由前文分析可知，求解方程组(1)时需求解湍流雷诺应力张量！可的湍流模型方程以及湍流热扩散系数&的方程％因此，四因子模型曰将分为两部分来阐述，即低雷诺数修正的k-模型与k『—模型％
1)低雷诺数修正的k-模型
标准k-模型是一种较为常见的二方程湍流模型，适合模拟充分发展的湍流流动，对低雷诺数过渡流和近壁区域计算结果不够理想％在四因子模型中，应用了针对低雷诺数修正的I-方程％与标准k-模型相似，引入了Bouss-inesq假设，认为湍流雷诺应力与牛顿剪切引力相似，均与应变呈正比，雷诺应力张量可写为湍流黏度与速度梯度的乘积，即：
YY=—空=—%("!+!T)M2k I(4)
%3式中，6为湍流雷诺应力％
与标准k-方程相似，低雷诺数修正的k-方程也定义了2个新的变量，即湍动能k及其耗散率)，如式(5)所示％
-=U'Z U')
与标准k-模型不同的是，本模型定义的湍流黏度为：
/=C$k G”(6)式中，C”为模型常数％为更加准确地模拟流体的近壁面行为，该模型定义了局部动力学特征时间t J91，由4个参数计算得到，即：
6”=f1A s+九")式(7)中4个参数均模化为归一化壁面距离R s与湍流雷诺数R et的函数，如式(8)所示％'f1$=(1—exp(—0.0714R a))z
A$=k&
,$(8)
f*=f$exp(—2.5X105Re t)
A$=5A1fl/Ret/i
归一化壁面距离R s和湍流雷诺数R e’定义如下：
\R S=!(/)0'25//
,(9)
[Re t=k t/-
式中：！为各网格单元到壁面的距离#为层流运动黏性系数％
为计算局部湍流黏度需得到式(7)中定 义的局部动力学特征时间t”，即需要解出归一化壁面距离和湍流雷诺数，进而需解出湍动能k及其耗散率)。修正后的k-方程如式(10)所示％
+▽・(U k)
▽•(/+比)k M P t—-
&+▽•(U s)=▽•(/+比)▽£+
(t'"J
C1-kP T-C-—f s
10)式中：
-P—=/(U+▽U T):V U
i f-=(1—exp(—0.3226R a))z•(11)
.(1—0.3exp(—0.0237Re())
1010原子能科学技术第55卷
=  1.5
C2e=  1.9
V(12)
C$=0.09
、"="k=  1.4
关于低雷诺数修正的I模型的更多阐释可参考文献)1*。
2)I-模型
如前文所述，通过式(2)将未知量湍流热流矢量［7于的求解转化为对湍流热扩散系数乂的求解。针对此，四因子模型定义了如下k-I-模型)2*。
与I模型类似，ki-模型同样定义了2个新变量：脉动温度方均值k-及其耗散率)-。
T-—2TY
V(13)
)-=/"T F
与SED模型不同的是，在四因子模型中将湍流热扩散率尙定义为式(14)所示的3个参数的乘积。
尙—C-kn-(14)式中：C-为模型常数；o-为局部热力学特征时间，用来更加精确地描述湍流流动中的换热特性「5*。
—f181-M+282-(15)式中，f i-、f2-、8-与8-4个参数被定义为归一化壁面距离R s、湍流雷诺数R et等变量的函数。
'f—(1—exp(—0.0526R S/槡&))・
(1—exp(—0.0714R J)
f282-—O u f2-r m c M
£2  1.3
诃叫槡槡&R愛4
81-=0.9k/)
(16)
-f2-—f-exp(—4U10—6@屛)(17)
I f2-—f-exp(—2.5X10-5R f)
式中,R为热力学特征时间与动力学特征时间的比值。
R—0—k(18)
O u k)-
综上所述，为计算湍流热扩散系数尙，需得到局部热力学特征时间G-，进一步需解出脉动温度方均值k-及其耗散率)-，总结kl-方程组如下：
(+▽・(U k-)=▽・(&M生)以-M
(\
T)2―)
V)+▽・(U-=▽・(&M生M
19)式中:C=1、C/i=0.925、C/2=0.9均为模型常数;R定义见式(11)/，与C2定义见式(20)。
C2=(1.9(1—0.3exp(—0.0237R#))—1)•V(1—exp(—0.1754R J)2
/-—T：£T)
(2)
2带绕丝棒束通道CFD模拟
2.1几何模型
本文选取德国卡尔斯鲁厄理工学院(KIT)盒间流实验「⑷中的带绕丝棒束组件作为CFD 模拟和模型验证对象。该棒束组件为六边形，包含有7根带绕丝燃料棒，每根棒表面施以均匀热流密度以模拟燃料棒对铅铋合金流体加热。从组件进口向出口看，7根绕丝顺时针旋转并缠绕在7根棒上。该带绕丝棒束组件具体几何模型与截面视图如图1所示，主要几何尺寸如下：/=20.5mm、"=16.0mm、？= 4.4mm、P=20.8mm、FF=61.05mm、L= 262mm、"h=10.34mm。
图1带绕丝棒束几何模型及其截面视图
Fig.1Geometric model and sectional view
ofwire-wrappedrodbundle
第6期 何少鹏等：基于OpenFOAM 的液态金属铅祕三维流动换热特性数值模拟研究1011
2. 2计算方法
本文选取绕丝螺距长度为262 mm 的1个 通道作为物理模型％进口温度均为573 K，施
加速度入口边界条件％出口为固定压力条件！ 组件壁和绕丝表面设为绝热壁面，棒束表面施
加均匀的热流密度％模拟中考虑了液态铅铋物
性随温度的变化，梯度项离散采用Gauss  linear  格式，扩散项采用Bounded  Gauss  linear 格式， 其余使用一阶迎风差分格式，分别使用SED 模型
与四因子模型进行计算，具体模拟工况列于表1% 2.3网格无关性分析
如图1所示，本文选取组件出口截面上1条 垂直线段A ，以沿线段A 的冷却剂温度变化特性
及速度分布为对象获得网格独立性解，如图2所
62840628
33222110T S .U I )^
网格无关性分析
Fig.2 Gridindependenceanalysis
2.4计算结果分析
带绕丝棒束通道内各截面处的速度、温度
云图如图3所示％结果表明，由于绕丝对流体
59077K
1-610
|L/|/(m-s _1)■-0.370.300.250.200.150.100.050.00
-600580570
图3各截面速度与温度云图
Fig.3 Velocityandtemperature
contourofseveralsections
示％从计算结果可看出，当网格数为166万、
196万和234万时，沿线段E 的冷却剂速度、温
度变化相差不大，故最终选取166万网格数的 网格为网格独立解开展后续数值模拟研究％
Table  1 Simulated  condition
表1模拟工况
工况
质量流量/(kg  • sT  )
入口速度/
(m  • s —1)P e(573 K )
R e(573 K )
10.716
0.0576
2 904
2
1.4320. 101525 8093
2.148
0. 15
2288 713
4
2.8640.20305
11 618
5  3.580
0.2538114522
注：加热功率为171 944 W/m 2
的搅混作用，各截面处速度分布呈现非对称、不 均匀特征，且沿轴向发生剧烈变化％由温度云 图还可看出，通道内不对称速度分布引起了不
对称温度分布。通道中心截面横流速度场矢量 图如图4所示，由于绕丝的搅混作用，横流速度
场呈现出旋转特点，且横流旋转方向与绕丝旋 转方向一致％
图4中间截面(z= 131 mm)横流速度矢量图
Fig.4 Cross  flow  velocity  vector  graph
ofmiddlesection  (K =131 mm
)