TOPAZⅡ反应堆本体流固共轭传热数值模拟

邹佳讯，郭春秋，孙 征

(中国原子能科学研究院 反应堆工程技术研究部，北京 102413)

TOPAZ-Ⅱ是由俄罗斯开发的热离子型空间核动力反应堆的堆型[1,2]，TOPAZ-Ⅱ反应堆不同于常规的陆上反应堆，其堆型要求相对紧凑、小型或微型化。TOPAZ-Ⅱ堆内的冷却剂为液态NaK合金，NaK共晶合金具有熔点低，热导率高等优点。 国内外针对TOPAZ-Ⅱ开发了一系列专用系统分析程序，如俄罗斯开发了ENSY系统分析程序对启动期间丧失冷却剂事故进行了分析，美国开发了热离子集成瞬态系统分析程序有TITAM，CENTRAR等[3～4]，国内开发的见诸报道的有西安交通大学开发的瞬态空间热离子反应堆系统分析程序TASTIN等等[5～6]，这些系统分析程序可以分析稳态或者瞬态工况下系统的整体参数及动态响应情况，关于TOPAZ-Ⅱ详细的堆本体三维流场温场方面的信息较少。堆芯入口流量分配数据是反应堆热工水力必不可少的计算输入，从反应堆力学分析的角度来看，堆芯燃料和部件可能会因为自身的温度梯度而产生热应力，流动和传热计算的到的温场及流场可为力学分析提供设计输入，一般情况下，固体的热应力分析和温度计算分析可能是耦合的，但由于应力形变相对于几何尺寸有时候是可以忽略的，所以在确定热应力之前可单独进行温度分布计算。目前的商用CFD软件无论是在国外还是国内，已越来越多应用于诸如水冷堆、液态金属冷却反应堆等反应堆热工水力行为的数值计算和研究中[7～9]。本文利用计算流体力学技术，针对TOPAZ Ⅱ反应堆本体进行流固共轭传热模拟，获得关键的热工水力参数及堆本体温度，为后续设计提供参考。

1 计算模型

1.1 TOPAZ-Ⅱ系统描述

该系统[10]由堆本体、电磁泵、辐射屏蔽(由钢重屏蔽及氢化锂轻屏蔽组成)、热排放系统、稳压系统以及重要的仪表控制、管道阀门等组成，如图1所示。其堆本体活性区内有氢化锆慢化剂，慢化剂上下为端部铍反射层，慢化剂外层为容器筒体；筒体外侧为侧铍反射层，侧铍反射层内有12个控制鼓其中3个为安全鼓，另外9个为调节转鼓，鼓体材料为铍，吸收体材料为碳化硼；活性区有 37个热离子燃料元件可以产生共计约4.5～5.5 kW电源，37个热离子元件分别坐落于均布在的慢化剂的37个竖直孔道内，孔道内由不锈钢内外套管组成单独的环形窄封流体通道，液态NaK合金经过37个环形窄封孔道竖直流动，不锈钢内外套管之间的间隙为冷却剂孔道，其中不锈钢内套管外径为24.5 mm,不锈钢外套管内径为25.9 mm，内外套管之间的间隙仅为0.7 mm。

37个流体通道两端为上下集流腔室，其中下集流腔上栅板、下集流腔下栅板与反应堆筒体焊在一起后形成冷却剂下集流腔；上集流腔上栅板、上集流腔下栅板焊在一起后形成冷却剂上集流腔。在慢化剂孔道内的冷却剂内外套管之间形成的流道与冷却剂上下集流腔相通，其中充满NaK冷却剂。图1(b)中1至37为燃料元件冷却剂通道所在位置编号。系统主要设计参数[1-6]如表1所示。

图1 TOPAZ-Ⅱ反应堆本体(停堆状态)Fig.1 Diagram of TOPAZ-Ⅱ reactor complex (shutdown condition)

表1 TOPAZ-Ⅱ的主要设计参数Table 1 Main design parameters of TOPAZ-Ⅱ

1.2 网格化分

以图2所示的TOPAZ-Ⅱ反应堆转换器为对象，可以采用分区的方法进行网格化法，上下上腔室分别为一个区域，中间37个冷却通道为第三个区域，中间的冷却通道区域两段分别用交界面与上腔室下栅格板、下腔室上栅格板相连。冷却剂上下腔室区域生成7层附面层贴壁网格。冷却剂周围的主要固体区域之间在区域设置上相互独立，以便灵活地的进行网格控制，固体区域在相贴的壁面上用交界面进行连接，用于数值模拟计算时热传导传递数据。计算结果与网格数目的敏感性方面，共划分了四套网格(见表2)进行计算，选择流体域NaK合金的最高温度和固体域慢化剂块的最高温度进行比较，从表中可以看出，网格总数在300万～400万的三套网格计算得到NaK合金最高温度、慢化剂最高温度之间的差距不到1 ℃，满足网格独立性的要求，后续给出的分析均为基于第四套网格的计算结果，第四套网格中具体的网格数目组成包括：流体域NaK合金2 802 616单元数；二氧化碳修补气体区域319 088单元数；慢化剂及端部铍反射层区域361 767单元数；侧铍反射层区域301 203单元数；碳化硼吸收体区域54 096单元数；控制鼓体区域180 516单元数。

图2 计算域网格示意图Fig.2 Mesh display of calculation zone

表2 数值结果网格独立性Table 2 Numerical results independence of numbers of mesh

1.3 材料物性

模拟计算中主要涉及的流体为NaK合金，固体有氢化锆慢化剂、铍反射层、不锈钢等，稳态数值模拟计算中需要的关键热物性热中，三者的密度分别为5 615 kg/m3、1 830 kg/m3和7 900 kg/m3，热导率分别为22.6 W/(m·K)、94.3 W/(m·K)、22.2W/(m·K)。流体换热中需要用到NaK合金的密度、比热容、热导率、动力黏度等热物性，固体的密度、比热容和热导率。冷却剂热物性参数(来自NaK工程手册)包括密度、热导率、比定压热容、动力黏度(见表3)，可以采取随温度的变化表进行插值或者用拟合函数载入求解。

表3 冷却剂热物性Table 3 Thermal and Transport Properties of Coolant NaK

碱性液态金属的对流换热数值模拟需要注意的是液态金属非常低的普朗特数Pr。

(1)

从表2可以得到NaK的Pr随着温度的变化曲线及拟合函数如图3所示，在高温区750～950 K范围内，Pr在0.006左右，非常低的Pr意味着分子热传导在整个传热过程中占据着较大的比重。

图3 NaK合金普朗特数随温度的变化拟合曲线Fig.3 Fitting curve of Pr of liquidNaK with temperature

CFX前处理器提供了多种湍流模型，如常见的标准k-ε模型、k-ω、SSTk-ω等雷诺平均纳维斯托克斯涡粘模型(RANS)，CFX前处理器中用自动近壁面处理方法(即壁面函数)对湍流流动中的近壁面的流动进行预测，而不需要精细化贴壁网格，在无量纲距离y+<300情况下，壁面函数均有效，且对y+没有最小值要求。表4为额定流量下部分数值结果针对湍流模型的敏感性，从中可以看出，三种常见的湍流模型得到的温度值相互之间的差距仅为0.2，但在反应流动和换热的特性数据来看，SST总体上更为接近经验关系式得到的值(见后文表5、表6)。

表4 部分数值结果对湍流模型的敏感性Table 4 Sensitivity of part results with certain turbulent models

后文CFD分析使用的是基于K-Omega的 SST模型，模型中的湍流普朗特数Prt(涡扩散系数和涡热传导系数的比值)默认值为0.85(FLUENT)/0.9(CFX)，适用于轻水或空气的模拟，不完全适用于模拟液态金属换热，有若干文献中提出或报道了计算Prt的经验关系式或推荐值，如下式：

(2)

Pe=Re×Pr

Prt=4.12，Pe<1 000(Xu Cheng[11])

(3)

对于本文中的模拟对象TOPAZ-Ⅱ反应堆液态NaK合金正常运行工况其雷诺数Re>4 000，其流动属于环管内的湍流流动，且其Pe数约30。利用上述Reynolds式计算得到的Prt约为4.1和Xu Cheng关系式中的4.12基本一致，因此需要更改软件的默认值0.9为4.1。

1.4 计算条件

上腔室入口采取质量流量入口边界，正常稳态额定工况总入口流量为1.3 kg/s，下腔室出口采取压力出口边界，设定为165000Pa，CFX中流体域中开启K-Omega SST湍流模型。物性数据用CFX的表达语言载入相应模块。活性区375 mm高度上的环形管道内壁面使用均匀热流密度，慢化剂释热率在径向上分为5个区域，每个区域的轴向归一化释热率分布如下图4所示，其中纵坐标为归一化释热率(针对慢化剂区域最大释热率进行归一化，慢化剂最大体积释热率约1.0×106W/m3)，分别对各区域的轴向分布进行曲线拟合，得到与轴向位置有关的释热分布，然后在前处理中使用表达式语言载入，而端部铍反射层(约0.1 kW)、侧铍反射层(约0.2 kW)、碳化硼(约0.4 kW)等其他固体区域发热相对于慢化剂(约4 kW)来说比较小，因此使用各区域平均体积释热率作为输入。侧反射层外表面设置辐射换热，发射系数为0.7，环境温度设置为300 K。

图4 慢化剂释热率分布Fig.4 Normalized heat generation rate distribution in moderator

3 计算结果与分析

3.1 流量分配

数值模拟计算完毕后，从CFD-POST提取37个燃料元件冷却通道的流量数据，通道冷却剂入口的流量分配系数由下式计算：

(4)

其中，εi、Mi分别为i通道的流量分配系数和流量。数值模拟得到流量分配结果如下图5所示。图中横坐标上的数字1～37代表各通道编号，编号与图1(b)中的号码一一对应，纵坐标为流量分配系数，由式(4)计算得到。尽管冷却剂从入口管进入上腔室这段流程流场较复杂，但由于37个燃料冷却通道首先入口流通面积相等，各入口面基本处于等压面，因此经分配后，各冷却剂管道内的冷却剂达到一相对均匀的状态。从图中可以看出，各通道流量分配因子在0.99～1.01；从堆芯中心至堆芯外围区域入口流量分配系数逐渐变大，各圈流量分配系数分布较为均匀。εi>1.0的燃料元件主要分布在堆芯外围区域，出口附近的几个通道(编号为20、21、36、37；30、31、32、33；24、25、26、27)因子较大。

图5 通道流量分配因子Fig.5 Mass-flow rate distribution factors

3.2 摩擦因子

为了验证流动数值计算的可靠性，从摩擦因子f的角度进行分析，摩擦因子f反映压降和流量的关系，其表达式如下：

(5)

其中ΔP/L为流动方向长度L上的压降，De为等效水力直径(通道特征尺寸)，ρ为流体密度，V为流体速度。对于光滑圆管湍流流动，经验关系式如下，其中Re为雷诺数：

(6)

德国学者V.Gnielinski[12]在光滑管湍流流动摩擦因子关系式基础针对环管湍流摩擦因子进行了修订，关系式如下：

fg=(1.8*lg10(Re*)-1.5)-2

(7)

(8)

(9)

式(8)即Gnienlinski增加的修订部分，α为环管内径di和外径do的比值，

为了对数值模拟得到的环管摩擦因子进行比较，在CFD-POST中的利用式(5)求得不同流量算例下环管的摩擦因子f_CFD，分别和利用经验关系式(6)和式(8)得到曲线进行比较，如图6所示，图中横坐标为额定流量65%～150%范围内雷诺数Re，纵坐标表对应的摩擦系数。从图中可以看出，数值模拟得到的摩擦系数与关系式(8)得到的摩擦系数fg之间的整体偏差不到5%，验证了流动计算的可靠性。

图6 不同雷诺数下摩擦系数比较Fig.6 Comparision of frcition factors under different Re number

3.3 努赛尔数

Nu是用于描述换热的一个无量纲数[式(10)]，针对圆管/环管内湍流流动换热，从二十世纪六七十年代至今，已有学者陆续提出计算液态金属强迫对流换热Nu若干经验关系式， 其一般表达形式见式(11)或其修订形式，在充分发展的流动换热流域中Nu一般随着Pe数的增加而增加，在换热稳定段约为定值。

(10)

Nu=Nuo+a×Pe

(11)

式中，a,b分别为经验关系式的拟合系数。

在CFD-POST中利用式(10)求得发热段入口一段距离上Nu的分布如图7所示，从图中看出其换热稳定段的Nu约为6.1。

图7 额定工况下Nu沿受热段轴向的分布Fig.7 Axial distribution of Nu number in inlet regions of NAK channel under nominal condition

各流量台阶算例下的努赛尔数的对比见表5。从表中可以看出，三者之间的差别很小，表明了CFX数值模拟换热分析的可靠性。

表5 各流量下努赛尔数Table 5 Nu under different mass flow rate conditions

3.4 温度场信息

环形通道内冷却剂温度轴向分布、慢化剂、端部铍反射层径向沿某点沿着轴向上的温度分布如图8所，图中纵坐标为归一化温度(局部温度/最高温度)，横坐标为归一化轴向位置(z/H)。冷却剂在活性区高度375 mm距离上温度成线性分布，在活性区外的两端为温度平台，分别于出入口温度接近，入口段温度为743 K的设置值，出口温度最高约846 K，仅与设计值843 K相差3 K。最高温度出现在慢化剂区域，在中心燃料通道附近区域，最高温度约为868 K，详细的剖面温度等值线见图9，10。

图8 慢化剂、端部铍反射层以及冷却剂轴向温度分布Fig.8 Axial temperature distribution of moderator、end reflectors and coolant

图9 横截面温度示意图Fig.9 Temperature distribution at axial cross section

图10 纵剖面温度示意图Fig.10 Temperature contour in the longitudinal cross section

4 结论

本文利用数值模拟方法对TOPAZ-Ⅱ反应堆进行了流固共轭传热计算分析，通过计算分别对反应流动和压降的摩擦因子以及反应传热的努塞尔数进行了比较，发现数值模拟计算得到的相应结果与使用经验关系式得到的计算结果非常接近，证明了数值模拟计算在反应堆流动与换热分析上的可靠性与正确性，所得到的详细流量分配数据可以为后续系统分析提供输入，部件温度三维场信息可以为力学计算提供接口，TOPAZ-Ⅱ的流固共轭传热数值模拟可以为后续相关堆型的优化设计奠定基础。