熔盐冷却球床堆热通道热工水力特性数值分析

牛 强 宋士雄 魏 泉 刘亚芬 梅龙伟 郭 威 蔡翔舟

熔盐冷却球床堆热通道热工水力特性数值分析

牛 强1,2宋士雄1,2魏 泉1,2刘亚芬1,3梅龙伟1,3郭 威1,3蔡翔舟1,3

1（中国科学院上海应用物理研究所 嘉定园区 上海 201800）2（中国科学院大学 北京 100049）
3（中国科学院核辐射与核能技术重点实验室 上海 201800）

基于计算流体力学(Computational Fluid Dynamics, CFD)通用计算程序Fluent，研究了模块化熔盐冷却球床堆(Pebble Bed Advanced High Temperature Reactor, PB-AHTR)中心热通道稳态热工水力行为。利用已开发的多孔介质流固两相局域非热平衡模型计算了球床堆中的压降、冷却剂的温场分布以及固相球床的温场分布，计算并比较了不同的多孔介质阻力因子(Ergun与KTA)对通道内的冷却剂流动以及温场分布的影响，并对丧失部分冷却剂情况下通道内的冷却剂及燃料温度进行了计算分析。结果表明使用不同的阻力因子对堆芯压降计算结果和流场的分布影响较大；而冷却剂温场及固相球床温场和球心的温度分布在不同的阻力因子下的差别较小，在PB-AHTR的设计参数下堆芯产生的热量能够被有效的输出，设计具有较大的安全裕度。计算结果对于球床堆的优化设计提供了一定的参考价值。

熔盐冷却球床堆(Pebble Bed Advanced High Temperature Reactor, PB-AHTR)，热工水力，多孔介质

2001年成立的第四代核反应堆国际论坛确立了包括超高温气冷堆、铅冷快堆、钠冷快堆、气冷快堆、超临界水堆和熔盐堆在内的第四代先进核反应堆类型，它们具有经济性好、防核扩散能力强、非能动的本征安全特性突出等特点，引起了核能界的广泛关注。美国UC Berkeley大学设计的900-MW模块化熔盐冷却球床堆(Pebble Bed Advanced High Temperature Reactor, PB-AHTR)采用包覆颗粒技术的燃料球作为燃料元件，使用7Li2BeF4熔融盐作为堆芯冷却剂，集中了第四代核能系统中熔盐堆和高温气冷堆各自的优势，具有高功率密度、低运行压力及热电效率高等特点，受到人们的普遍青睐[1-3]。

UC Berkeley大学设计的PB-AHTR方案中的一种堆芯设计参数[4]为：活性区高度3.2 m，包括7个边心距62.5 cm的六棱柱状组件，每个组件有19个直径19.8 cm的燃料球填充孔道，此孔道同时也作为冷却剂通道。孔道长度为3.2 m，填充3 286个直径dp为3 cm的燃料球，每个截面上分布31个燃料球，平均填充因子（孔隙率）ε=47.14%。每个孔道冷却剂的质量流率6.226 7 kg·s-1，入口温度873 K，堆芯出口平均温度为1 173 K。图1、2为堆芯结构及冷却剂通道内燃料球的填充示意图[4]。

图1 堆芯模型Fig.1 Reactor core model.

图2 通道内燃料填充图Fig.2 Fuel pebble packing model.

本文基于CFD通用计算程序Fluent[5-6]，利用已经开发的局域非热平衡多孔介质模型[7-10]，对PB-AHTR堆中心热通道内的热工水力特性进行了分析，研究了不同情况下通道内冷却剂的流动行为、冷却剂的温场分布以及固相球床的温场分布。

1 理论模型及计算方法

由于高温球床堆堆芯热通道具有强内热源，冷却剂与固相球床之间的温度梯度较大，Fluent自带的基于局域热平衡假设的多孔介质模型在这种工况下不再适用[8]。本文通过在基于表观速度的动量方程上附加阻力源项Si来表征多孔介质区域的流动行为。冷却剂的温场通过流体能量场方程求解，固相球床的温场利用Fluent的UDS (User Defined Scalar)功能求解，在迭代过程中通过交换流体能量场方程与UDS的源项实现两相的热耦合。通过与稳态基准题比较，显示本模型对球床堆堆芯稳态主要热工参数计算结果与国际通用的球床堆热工程序THERMIX、TINTE符合较好(表1)，可以用于球床式反应堆的热工水力分析[10]。

相对于传统的球床堆热工分析程序，CFD方法对计算网格的依赖程度较高，为了得到合理的计算结果，必须进行网格敏感性分析。在相同的边界条件与求解精度下，对堆内中心热通道建立了三套网格进行计算。表2为主要热工水力参数（流体温度、固相温度等）在不同的网格密度下的计算值。在网格体密度为449 445时各主要计算参数与最密的网格(717 020)下的结果无明显变化(＜5%)，已达到网格无关性的要求。

表1 Fluent结果和THERMIX、TINTE的比较[10]Table 1 Calculation result comparison between CFD and THERMIX, TINTE[10].

表2 网格无关性分析Table 2 Grid independent analysis.

在Fluent中采用SIMPLE算法求解局域非热平衡多孔介质模型的流体控制方程，模型的控制方程如下：

连续性方程(,,,=i jx y z)：

动量守恒方程：

式中，Si为Ergun[11]阻力因子或德国KTA[12-13]推荐的应用于高温气冷堆计算的阻力因子，其中Ergun阻力因子对于球床式反应器具有广泛的适用性，多用于固定床和流化床反应器流体分析计算，KTA推荐的阻力因子与高温气冷堆实验值符合较好。对于熔盐冷却球床堆，目前尚无实验验证其阻力因子的适用性，本文考虑了这两种不用的阻力因子对热工计算的影响，在缺乏实验数据的初步设计时可提供有用的参考。两种阻力因子的表达式分别在式(3)、(4)中给出：

式中，μ为熔盐的动力学粘度；ρ为流体密度；u为表观速度；dp为燃料球的直径；ε为孔隙率。

流体温场和固相球床温场分别通过能量守恒方程(5)、(6)进行求解：

其中，流固两相之间的换热量Ss为：

对流换热系数可以用球床传热努赛尔数(Nu)与球床几何参数求出：

Nu有两个常用的经验关系式，分别为Wakao[14]式(7)与KTA推荐的应用于高温气冷堆的关系式(8)，其中Wakao公式多用于化学工程反应器的热工分析，KTA推荐的Nu数计算公式是通过球床式高温气冷堆的实验数据拟合而成的。

其中，z为轴向高度，单位：cm。

考虑到通道径向孔隙率分布的不均匀性，本文中距孔道中心以7.5 cm为半径的圆内设置孔隙率ε=0.395 4、圆外ε=0.573 8。

7Li2BeF4熔盐物性和热工参数[15-16]如下：

球床的热传导系数的计算采用Zehner-Schlunder[17]模型。

流体的雷诺数Re约4 500，在多孔介质区域为完全的湍流，计算时开启标准k-ε湍流模式。在迭代过程中交换源项Ss=-St实现温场的耦合计算。

2 堆芯热通道热工水力特性分析

PB-AHTR热通道区域高度为H=3.2 m，操作压力为大气压，入口温度为873.15 K，质量流率m˙=6.226 7 kg·s-1，折算入口速度为0.101 76 m·s-1。本节冷却剂与球床的对流换热系数采用Wakao关系式，分别对附加Ergun阻力源和KTA阻力源下的堆芯热通道内热工水力特性进行分析。

2.1 压降分析

图3、4分别为Fluent附加Ergun阻力源和附加KTA阻力源得到的堆芯压降分布。堆芯在附加Ergun阻力源下的压降为19 370.2 Pa，附加KTA阻力源计算得到的压降为15 199.3 Pa，与Ergun阻力源得到的压降有较大差别，使用不同的阻力源经验公式对堆芯压降计算结果的影响较大。堆芯的压降关系到反应堆的经济性能评估与各类工况下回路的唧送功率匹配，因此在进行熔盐冷却高温堆堆芯设计时，需要借助相关热工水力学实验得到适合于特定堆型的阻力源经验公式。

图3 Ergun阻力因子压降Fig.3 Pressure drop with Ergun factor.

图4 KTA阻力因子压降Fig.4 Pressure drop with KTA factor.

2.2 流场分析

附加Ergun阻力源下冷却剂在Z平面垂直方向上速度场分布如图5、6所示，根据计算可得通道中心位置处流速最低降到0.075 m·s-1，这是因为靠近通道中心区域的孔隙率要高于外环区域所致，对流速影响最大的是通道内的孔隙率分布。

图5 Ergun阻力源项X=0平面Z方向流速Fig.5 Z-velocity at X=0 plane with Ergun factor.

图6 Ergun阻力源项Z平面Z方向流速Fig.6 Z-velocity at Z-plane with Ergun factor.

由图7、8可见，KTA阻力因子计算的Z方向流速的最小值要低于利用Ergun因子计算的流速，流场在1.8 m以下的区域分布较为均匀。由于多孔介质的阻力源项较大，不论是附加哪种阻力因子，靠近壁面处的流速只有在入口处的壁面效应比较明显，壁面处的流速分布在Z轴方向上也不是很均匀。本文采用的宏观多孔介质模型，在流动方向上存在较大的阻力源，其对流速的影响远大于边界层的影响，因此速度场在近壁面处无明显壁面效应。

图7 KTA阻力源项X=0平面Z方向流场Fig.7 Z-velocity at X=0 plane with KTA factor.

图8 KTA阻力源项Z平面Z方向流场Fig.8 Z-velocity at Z-plane with KTA factor.

2.3 冷却剂和固相球床的温场分析

根据程序计算，在附加Ergun阻力源情况下的中心热通道总功率为5 396 041 W，通道内熔盐平均温度为1 057.37 K，最高温度为1 266.26 K。固相球床平均温度为1 152.46 K，最高温度为1 348.65 K。球床与冷却剂的平均温差在100 K左右。熔盐出口平均温度为1 230.42 K。图9为通道内X=0与Y=0截面上冷却剂的温度分布，图10为X=0与Y=0平面上固相球床的温度分布。图11为堆芯Z=0.4 m、1.6 m和2.6 m处固相球床温场分布，中心区域与边缘区域温差达到了80 K以上。

图9 冷却剂温场分布Fig.9 Coolant temperature profiles.

图10 球床温场分布Fig.10 Pebble bed temperature profiles.

图11 Z平面球床温度分布Fig.11 Pebble bed temperature at Z-plane.

附加KTA阻力源项下的通道内熔盐平均温度为1 057.74 K，最高温度为1 262.88 K。固相球床平均温度为1 153.67 K，最高温度为1 348.68 K。熔盐出口平均温度为1 230.93 K。图12为堆芯通道内X=0与Y=0平面上冷却剂温度分布，图13为X=0与Y=0平面上固相球床温度分布。图14为堆芯Z=0.4 m、1.6 m和2.6 m处固相温场分布，中心区域与边缘区域温差达到了60 K左右。

由图9和12可知，Ergun和KTA阻力源对于冷却剂的温场分布的影响差别较小。图10和13表明在堆芯中心热通道的温度梯度较大，这是因为规则排布的球床由于受到通道本身尺寸的限制，使孔隙率分布在径向有较大的不均匀性，从而造成了较大的温差。通道内孔隙率分布对于通道内冷却剂及固相球床的温场分布有很大影响，在装载燃料球时应尽量使通道壁面处的孔隙率低于通道中心，有利于冷却通道的石墨壁面，延长其使用寿命。

图12 冷却剂温场分布Fig.12 Coolant temperature profiles.

图13 球床温场分布Fig.13 Pebble bed temperature profiles.

图14 Z平面球床温度分布Fig.14 Pebble bed temperature at Z-plane.

仿真计算表明，使用不同的阻力因子对堆芯压降计算结果影响较大，对冷却剂温场及固相球床温场分布影响较小。在堆芯初步设计时，由于缺乏实验数据，阻力因子的选取具有一定的不准确性，但是在稳态工况下，冷却剂与球床的温场分布与阻力因子的相关性不大。通过比较两种阻力因子对堆芯主要热工参数的影响可认为，在进行初步设计的热工分析时，可以选取其中一种阻力因子而不影响主要分析结果。对于堆芯经济型的评估与回路唧送功率计算时，阻力因子的选取具有很大影响，未来可通过冷态的堆芯流体实验验证阻力因子的适用性。

3 失去部分冷却剂情况下热通道内热工水力行为分析

球床堆一个很重要的安全问题是堆芯在各种工况下的燃料球温度是否超过了安全阈值。在稳态运行工况下，燃料球的温度不能超过1 250 °C，在失冷事故工况下燃料球的温度不能超过1 600 °C[15]。本节分析了在部分失冷工况下PB-AHTR堆芯热通道的冷却剂及球床的温度分布，为优化堆芯设计及安全分析提供参考数据。由于中心单通道的热工水力-中子学反馈对通道的功率分布及功率水平影响较小，因此在稳态分析时未考虑物理热工反馈效应。

通过调节热通道的入口速度来调节冷却剂流量，计算失去10%(case 1)、20%(case 2)、30%(case 3)、40%(case 4)、50%(case 5)冷却剂情况下通道内的冷却剂以及固相球床的温度分布。阻力源项采用KTA推荐的关系式，并比较了在KTA与Wakao关系式下热工参数的区别。图15、16显示了堆芯熔盐、球床的温度随不同冷却剂流量的变化关系。

图15 利用Wakao公式计算部分失冷工况温度Fig.15 Temperature of loss of coolant with Wakao equation.

图16 利用KTA公式计算部分失冷工况温度Fig.16 Temperature of loss of coolant with KTA correlation.

图15 、16显示随着冷却剂流量的丧失，堆芯热通道内的熔盐及固相球床的温度上升很快。但是在失去50%冷却剂的情况下，球床的最高温度仍然没有超过燃料球的限制温度。在PB-AHTR的堆芯结构设计中冷却剂通道上部为一个装料的球床区域，通道的冷却剂分配也通过这个区域实现。本文计算结果表明，只要堆芯热通道的冷却剂入口区域设计合理，不产生较大区域的阻塞，则堆芯的燃料温度就可以保持在安全阈值之下。从图15、16中还可以发现，利用不同的换热系数计算得到的热工参数基本没有差别。这与文中PB-AHTR的冷却剂流速较低有关，在较低的雷诺数下，Wakao公式与KTA关系式的计算值相差不大。结果表明，在低流速稳态工况下利用宏观的多孔介质模型计算的温场对于换热系数的选取不是很敏感。

4 结语

本文基于CFD通用计算程序Fluent，利用开发的局域非热平衡的多孔介质模型，研究了模块化PB-AHTR反应堆中心热通道内热工水力特性，计算了冷却剂和固相球床的温场分布，分析了堆芯流场及特定轴向高度上的温场分布。计算并比较不同的多孔介质阻力因子(Ergun与KTA)下通道内的流动以及温场分布。结果表明：(1) 采用不同的阻力因子对通道的压降与流场分布影响较大。(2) 通过对附加Ergun以及KTA阻力因子的计算结果分析表明，冷却剂平均温度、最高温度以及固相球床的平均温度、最高温度在不同的阻力因子下差别均很小。KTA推荐的阻力因子适用范围0.36＜ε＜0.42，其适用性有待实验检验，Ergun阻力因子的适用性较强，但是对于新型熔盐冷却球床堆的适用性也有待实验检验。Ergun与KTA阻力因子都有可调节的常数因子，可利用实验数据拟合合适的常数因子，用于计算预测通道内的热工水力特性。(3) 最后分析了堆芯热通道在丧失部分冷却剂情况下的热工水力学特性，确认了在PB-AHTR的设计参数下堆芯产生的热量能够被有效地输出，并且在失去部分冷却剂的工况下能够保证反应堆仍处于安全阈值之内，PB-AHTR的设计具有较大的安全裕度。

1 Ball S J, Forsberg C W. Advanced High-Temperature Reactor (AHTR) loss of forced circulation accidents[R]. USA: Oak Ridge National Laboratory, 2004

2 Bardet P, Blandford E, Fratoni M, et al. Design, analysis and development of the modular PB-AHTR[R]. USA: University of California, Berkeley, 2008

3 Fratoni M. Development and applications of methodologies for the neutronic design of the Pebble Bed Advanced High Temperature Reactor (PB-AHTR)[R]. USA: University of California, Berkeley, 2008

4 Orzáez J A. Neutronics analysis of a modified pebble bed advanced high temperature reactor[D]. USA: Ohio State University, 2009

5 Fluent 6.3 User’s Guide[OL]. Fluent Inc. http://aerojet. engr.ucdavis.edu/fluenthelp/html/ug/main_pre.htm, 2006

6 Fluent 6.3 UDF Manual[OL]. Fluent Inc. http://aerojet.engr.ucdavis.edu/fluenthelp/html/udf/node2. htm, 2006

7 闫晓, 肖泽军, 黄彦平, 等. 多孔介质中流动换热特性的研究进展[J]. 核动力工程, 2006, 27(1): 78-82

YAN Xiao, XIAO Zejun, HUANG Yanping, et al. Research progress on flow and heat transfer in porous media[J]. Nuclear Power Engineering, 2006, 27(1): 78-82

8 杜建华, 王补宣. 带内热源多孔介质中的受迫对流换热[J]. 工程热物理学报, 1999, 20(1): 69-73

DU Jianhua, WANG Buxuan. Research on the forced convective heat transfer for fluid flow through porous media with internal heat source[J]. Journal of Engineering Thermophysics, 1999, 20(1): 69-73

9 Roy C T, Simon T, David G, et al. An initial non-equilibrium porous-media model for CFD simulation of stirling regenerators[R]. USA: National Aeronautics and Space Administration, 2006

10 宋士雄, 魏泉, 蔡翔舟, 等. 基于CFD方法的球床式高温气冷堆稳态热工水力分析[J]. 核技术, 2013, 36(12): 120601

SONG Shixiong, WEI Quan, CAI Xiangzhou, et al. High temperature gas-cooled pebble bed reactor steady state thermal-hydraulics analyses based on CFD method[J]. Nuclear Techniques, 2013, 36(12): 120601

11 Ergun S. Fluid flow through packed columns[J]. Chemical Engineering Progress, 1952, 48(2): 89-94

12 KTA 3102.3. Reactor core design of high-temperature gas-cooled reactors Part 3: loss of pressure through friction in pebble bed cores[S]. 1981

13 KTA 3102.2. Reactor core design of high-temperature gas-cooled reactors Part 2: heat transfer in spherical fuel elements[S]. 1983

14 Wakao N, Kaguei S, Funazkri T. Effect of fluid dispersion coefficients on particle-to-fluid heat transfer coefficients in packed beds: correlation of nusselt numbers[J]. Chemical Engineering Science, 1979, 34(3): 325-336

15 de Zwaan S J. The liquid salt pebble bed reactor[D]. Netherlands: Delft University of Technology, 2005

16 David S. Molten salt coolants for high temperature reactors[R]. IAEA, 2009

17 Zehner P, Schlunder E U. Thermal conductivity of granular materials at moderate temperatures[J]. Chemie Ingenieur Technik, 1970, 42: 933-941

CLC TL349

Thermal-hydraulics numerical analyses of Pebble Bed Advanced High Temperature Reactor hot channel

NIU Qiang1,2SONG Shixiong1,2WEI Quan1,2LIU Yafen1,3MEI Longwei1,3GUO Wei1,3CAI Xiangzhou1,3
1(Shanghai Institute of Applied Physics, Chinese Academy of Sciences, Jiading Campus, Shanghai 201800, China)
2(University of Chinese Academy of Sciences, Beijing 100049, China)
3(Key Laboratory of Nuclear Radiation and Nuclear Energy Technology, Chinese Academy of Sciences, Shanghai 201800, China)

Background: The thermal hydraulics behavior of the Pebble Bed Advanced High Temperature Reactor (PB-AHTR) hot channel was studied. Purpose: We aim to analyze the thermal-hydraulics behavior of the PB-AHTR, such as pressure drop, temperature distribution of coolant and pebble bed as well as thermal removal capacity in the condition of loss of partial coolant. Methods: We used a modified FLUENT code which was coupled with a local non-equilibrium porous media model by introducing a User Defined Scalar (UDS) in the calculation domain of the reactor core and subjoining different resistance terms (Ergun and KTA) to calculate the temperature of coolant, solid phase of pebble bed and pebble center in the core. Results: Computational results showed that the resistance factor has great influence on pressure drop and velocity distribution, but less impact on the temperature of coolant, solid phase of pebble bed and pebble center. We also confirmed the heat removal capacity of the PB-AHTR in the condition of nominal and loss of partial coolant conditions. Conclusion: The numerical analyses results can provide a useful proposal to optimize the design of PB-AHTR.

Pebble Bed Advanced High Temperature Reactor (PB-AHTR), Thermal-hydraulics, Porous media

TL349

10.11889/j.0253-3219.2014.hjs.37.070602

中国科学院战略性先导科技专项资助项目(No.XDA02010200)资助

牛强，男，1987年出生，2011年毕业于洛阳师范学院，现为硕士研究生，研究方向为反应堆热工水力

蔡翔舟，E-mail: caixiangzhou@sinap.ac.cn

2014-03-05，

2014-04-16