基于OpenFOAM的乏池内流动与换热计算

Calculation of Flow and Heat Transfer in Spent Pool Based on OpenFOAM

摘要: 通过OpenFOAM对湿法贮存下的乏燃料池内的热工水力问题进行仿真分析，得到以下结论：在燃料组件区，由于多孔介质下的阻力，该部分流场流速较小。位于进水口正下方的套筒内的水流受到强迫对流的影响产生了自上而下的流动，远离入水管正下方的贮存格架内的水流受到自然对流的作用，均自下而上地流动。这种流动导致了经由注水直接冲刷的正下方的格架内的水温极低，邻近套筒内的局部高温；入水管正下方的格架内的温度分布受到水流方向的影响，呈现出底部温度高、顶部温度低的分布特征。由于位于出口下方的套筒与邻近套筒组成了回路，使得周围邻近套筒的水温出现极大值。

Abstract: Through OpenFOAM, the thermal and hydraulic problems in the spent fuel pool under wet storage are simulated and analyzed, and the following conclusions are drawn: In the fuel assembly area, due to the resistance in the porous medium, the flow rate of the shunt field is small. The water flow in the sleeve directly below the inlet flows from top to bottom under the influence of forced convection, and the water flow in the storage lattice far away from the inlet pipe directly below the natural convection flows from bottom to top. This flow results in extremely low water temperature in the lattice directly below the flush by water injection and local high temperature in the adjacent sleeve; The temperature distribution in the grid directly below the inlet pipe is affected by the direction of water flow, showing the distribution characteristics of high temperature at the bottom and low temperature at the top. As the sleeve located below the outlet forms a loop with the adjacent sleeve, the water temperature of the surrounding adjacent sleeve is maximized.

1. 引言

无论从装机规模还是核电机组的制造与运营，中国的核电技术已经达到国际先进水平 [1]。但是作为核电的产物，乏燃料在相当长的时间内具备放射性，必须对其进行谨慎处理 [2] [3]。湿法储存长期以来风险较小，是许多核电厂选择乏燃料中间储存的主要方式，然而日本福岛事故的发生对现有的湿法储存的安全性提出了质疑。当安全保障措施不够时，极易发生重大事故 [4]。

乏燃料贮存水池作为过渡暂时性的燃料存储，其设计、运行、管理直接影响整个核电厂的运行安全。对乏池内的热工水力问题(目标问题)进行三维的计算流体力学(CFD)仿真是十分必要的，其可提供详细的流动传热细节和对关注区域的可靠预测。

文献 [5] [6] 研究了模块式高温气冷堆的非能动余热排出系统，在紧急事故后，堆芯周围的水冷壁可以将乏燃料的热量带走。文献 [7] [8] 也研究了乏燃料池内的几种能动冷却方案，都是大同小异。比较有名的是美国基于AP1000非能动设计的SFS系统，在断电的极端情况下，非能动安全壳冷却水箱及乏燃料清洗池可为水池补水，依靠池水沸腾蒸发来短期冷却乏燃料水池 [9]。

上述都对乏燃料的冷却进行了多方面的研究，但是很少有在开源软件OpenFOAM的基础上进行的仿真设计，针对乏池内的自然对流的优化设计，很难做到前后处理以及求解器的二次开发，高效地进行仿真工作。

本文针对水作为冷却剂的、盒式组件构成的乏池的热工水力问题进行仿真，利用开源软件OpenFOAM进行仿真计算，主要关注冷却剂的单相湍流和换热问题，求解燃料棒的温度和冷却剂的温度、速度等参数，以实现对实际乏池设计和优化的有效指导和借鉴。本文研究了湿法贮存下的乏燃料池内的热工水力问题，研究对象由22台贮存格架与2台设备组成。

2. 理论与模型

2.1. 多孔介质模型

采用多孔介质模型的能量方程如下：

$\frac{\partial }{\partial t}\left(\gamma {\rho }_f{E}_f+\left(1-\gamma \right){\rho }_s{E}_s\right)+\nabla \cdot \left({\rho }_{f}u{H}_f\right)=\nabla \cdot \left(k_{\text{eff}}\nabla T+\tau \cdot u\right)+S_f^h$ (1)

这里 $E_f$ 为流体总内能(包括动能)， $H_f$ 是流体总焓值， $E_s$ 固体总内能， $k_{\text{eff}}$ 为有效热导率， $S_f^h$ 是源项。而有效热导率是流体和固体热导率的体积平均：

$k_{\text{eff}}=\gamma k_f+\left(1-\gamma \right)k_s$ (2)

2.2. 自然对流模型

对于稳态条件下的自然对流控制方程，对于工质为水的情况，乏池内温度变化较小，可以采用Boussinesq假定：

$\rho ={\rho }_{ref}\left(1-\beta \left(T-T_{ref}\right)\right)$ (3)

在该假设下，有如下关系式

$\nabla p-\rho g=\nabla p_{rgh}+gh\nabla \rho$ (4)

这里 $p_{rgh}=p-\rho gh$ ； $gh=g\cdot r$ 表征的是垂直高度与重力加速度的乘积。

连续性方程为

$\nabla \cdot \left(\rho U\right)=0$ (5)

动量方程为

$\nabla \cdot \left(\rho UU\right)-\nabla \cdot \left(\mu \nabla U\right)=-\nabla p_{rgh}+gh\nabla \rho$ (6)

在能量方程中，假定忽略流体的动能，并认为热通量占主导地位，忽略由应力和重力引起的相关能量交换，那么可以写为：

$\nabla \cdot \left(\rho UT\right)=\nabla \cdot \left({\alpha }_{\text{eff}}\nabla U\right)$ (7)

这里 ${\alpha }_{\text{eff}}$ 为有效热扩散率。在多孔介质区域就是工质与固体材料的平均热扩散率。

${\alpha }_{\text{eff}}=\gamma {\alpha }_f+\left(1-\gamma \right){\alpha }_s$ (8)

3. 几何与网格划分

计算区域为三维的多通道燃料组件，结构比较复杂，下图截取单通道部分，主要包括上管座、下管座、端部格架(2个)、大格架(6个)和小格架(3个)、导向管和燃料棒，如图1所示。

Figure 1. CAD diagram of geometric model

图1. 几何模型CAD图

根据乏燃料贮存池的特点，本文将乏燃料池的多通道燃料组件看做多个矩形体组装而成。背景网格分为上下两层，下层由于结构较为复杂，所以平面内的网格尺寸较小。流体区域为：燃料组件所在的多孔介质区域，燃料组件与筒壁间的狭窄内流道区域，燃料棒间的外流道区域，底板所在的多孔介质区域，和底座周围的流域。

网格采用代码进行自动划分，部分代码如下图2：

Figure 2. Code for grid setup

图2. 网格划分设置的代码

因为冷水补给造成的水面波动很小，所以可将乏燃料贮存池与空气接触面视为同一个水平面 [10]。在该平面区域附近的传热受自然蒸发和自然对流双重影响。为了计算更接近于真实环境，本文利用等效传热系数法来计算在自由水池表面处的传热。

乏燃料池进行冷却时，冷却系统涵盖了自然对流和强迫对流两种情况，本文主要研究乏池的自然对流的冷却情况。当前很多研究人员针对自然对流的特性进行了CFD分析。根据之前的研究发现：在乏池的自然对流中，需要考虑流场和温度场强烈耦合，因此在计算过程中经常收敛较困难。特别乏池的自然对流在层流向湍流过渡的过渡区域时，收敛非常困难 [11] [12]。所以，考虑到收敛问题，本文模型的网格都与流动方向基本一致，采用结构网格进行处理，减少计算的假扩散现象。

4. 求解

该求解器的运行环境是Linux，推荐使用Ubuntu 16.04。另外所要求的环境有OpenFOAM v5.0，GNU GSL v2.5与GNU M4 v1.4.17。

在前处理阶段，需要根据自身需求对乏池的结构进行建模与网格划分，然后需要进行相关OpenFOAM字典文件的编写。

在进行预估场的计算时，OpenFOAM中的potentialFoam可以对初始的速度场及压力场进行预估，从而减少后续求解时的计算压力，减少求解发散的风险。本文使用k-OmegaSST模型进行计算。

4.1. 初始条件和边界条件

初始条件与边界条件定义如下表。入口为速度入口，−z方向2 m/s，温度为293 K；出口为压力入口，p = 101325 Pa；顶部壁面为滑移壁面，温度为293 K；其他壁面为无滑移壁面，湍动能、湍流耗散率均采用壁面函数，具体设置如下表1。

Table 1. Boundary condition setting

表1. 边界条件设置

对于薄板边界，compressible::thermalBaffle1D< hConstSolidThermoPhysics >，具体定义内容如下图3所示，表2为薄板模型的属性。

Table 2. Solid attribute parameters of thin plate model

表2. 薄板模型固体属性参数

Figure 3. Definition of temperature boundary of thin plate

图3. 薄板的温度边界定义

refValue等都是与薄板换热相关的参数。对于1D薄板的某侧p的换热有，

$Q_r+\frac{k_w}{\delta }\left(T_i-T_p\right)=\frac{k_s}{{\delta }_s}\left(T_p-nbr{T}_p\right)$ (11)

其中，

$Q_r$ 是辐射换热；

$k_w$ 是靠近壁面的流体的热导率；

$\delta$ 是壁面面心到网格中心的距离；

$k_s$ 是薄板的热导率；

${\delta }_s$ 是薄板的厚度；

$T_i$ 是该侧附面层网格中心的流体温度；

$T_p$ 是该侧壁面温度；

$nbr{T}_p$ 是薄板另侧壁面温度。

使用Darcy-Forchheimer公式对多孔介质区流阻进行建模，

$S_i=-\left(D_{ij}+\frac{1}{2}\left|u_{kk}\right|F_{ij}\right)u_i$ (12)

D为一次项系数，F为二次项系数。该字典支持自定义坐标系，本算例中，指定了原点origin为(0, 0, 0)，e1 (1, 0, 0)，e2 (0, 1, 0)，即与原先坐标系保持一致。

fvSchemes字典中指定时间项差分格式为steadyState，扩散项采用中心差分Gauss linear，对流项均采用一阶迎风差分格式bounded Gauss upwind。

fvSolution字典中指定了线性求解器的参数、收敛判据和松弛因子。相对残差relTol均取得较大，以保证稳定性；优化参数nCellsInCoarsestLevel取平均网格数的0.5次方。收敛判据为1e⁻³。对U，h，k，epsilon取0.5的松弛因子；对p_rgh场取0.3；对rhok场取0.001。

4.2. 自定义热源分布

通过编写新的源项类scalarVerticalProfileSource，借助于批处理脚本，实现在csv表格文件中自定义热源分布的功能。功能包括默认均匀热源强度分布、用户指定均匀热源强度分布、用户指定竖直方向非均匀热源强度分布。

打开位于preProc/config下的heSourceProfile.csv文件，进行相关乏棒热源强度的指定，该文件样式如下：

第一行为默认均匀热源强度分布的指定；第二行为自定义指定的开始行，不会被脚本读取；第三行为自定义均匀热源强度分布的指定方式，按顺序填入乏棒的编号，以逗号分隔，最后填写所需的热源强度即可；第四行为自定义竖直方向非均匀热源强度分布的指定方式，在乏棒编号后按顺序填入从底部到顶部的热源强度等距分布，以空格分隔，数据个数不限，脚本将自动读取该数据系列生成字典文件，而后求解器将自动获取乏棒cellZone底部和顶部的坐标，并通过插值对相关cellZone内每一个网格的源项进行指定。对于未指定的乏棒，脚本会自动采用默认均匀热源强度分布。

脚本读取CSV文件，在case/system下的fvOptions字典文件中自动生成所有乏棒的热源分布，样式如图4。

通过在求解器中添加相关语句，实现了源项的可视化，此例的源项可视化如左图所示，最左侧乏棒有指定的均匀热源强度，中间5个乏棒有默认的均匀热源强度，最右侧乏棒有指定的正弦分布热源强度(如图5所示)。

5. 计算结果

图6给出的是一区某格架内薄板给流场带来的阻碍作用。可以看出在区间0.14~0.32是燃料组件区，由于多孔介质下的阻力，该部分流场流速较小；区间0.32~0.35给出了是套筒内流道的流速情况；区间0.35~0.4给出的是套筒与套筒间的外流道的流速情况。

Figure 4. Automatic generation of heat source distribution code

图4. 自动生成热源分布代码

Figure 5. Effect diagram of heat release rate distribution loading

图5. 释热率分布加载效果图

Figure 6. Velocity distribution of several groups of sleeves along the Y-axis

图6. 若干组套筒沿y轴方向的流速分布

以下分别给出在若干典型的平面内得到的流场与温度云图。

1) 平面z = 3 m，温度云图如下图7。

Figure 7. Temperature distribution in the z = 3 m plane

图7. z = 3 m平面的温度分布

从图中可以看出，位于进水口正下方的套筒内的水流受到强迫对流的影响产生了自上而下的流动，这种流动导致了经由注水直接冲刷的正下方的格架内的水温极低，邻近套筒内的局部高温；而其他区域的套筒的水流在自然对流的作用下自下而上的流动。

2) 包含速度入口的y = 0平面，其温度云图如图8。

入口截面内的温度云图表面，温度极大值与极小值之间的差约为10 K。且入水管正下方的格架内的温度分布受到水流方向的影响，呈现出底部温度高，顶部温度低的分布特征。而其他部位的水温呈现出自然对流条件下的分布特征，即顶部温度相对较高。

Figure 8. Temperature distribution in y = 0 plane

图8. y = 0平面的温度分布

从图8中可以看出，位于出口下方的套筒与邻近套筒组成了回路，使得周围邻近套筒的水温出现极大值。

从图9中看出，由入口进入乏池的水流自上而下的流动，并且位于入水管正下方的贮存格架内的水流受到强迫对流的影响，也是自上而下流动。而远离入水管正下方的贮存格架内的水流受到自然对流的作用，均自下而上的流动。而位于格架上方的乏池内，出现了大范围的涡旋结构。

Figure 9. Velocity distribution in y = 0 plane

图9. y = 0平面的速度分布

3) 远离速度入口的x = 1.16 m平面，其云图为图10~12。

从图10中可以看出，因为乏燃料的衰变热的不断释放，水温在乏燃料贮存格架密集处较高，同时由于远离了入口，使得套筒中的流体进行自然对流，温度自下而上逐渐升高。

从图11、图12中可以看出，速度梯度较大的区域其湍动能也较大，其换热也更加强烈，热量源源不断的被带走。在入口的大空间，由于强烈的换热，从而水温更加均匀。

Figure 10. Temperature distribution in x = 1.16 m plane

图10. x = 1.16 m平面的温度分布

Figure 11. Velocity in x = 1.16 m plane

图11. x = 1.16 m平面的速度

Figure 12. Turbulent kinetic energy in x = 1.16 m plane

图12. x = 1.16 m平面的湍动能

6. 结论

通过OpenFOAM对湿法贮存下的乏燃料池内的热工水力问题进行仿真分析，得到以下结论：

1) 在燃料组件区，由于多孔介质下的阻力，该部分流场流速较小。

2) 位于进水口正下方的套筒内的水流受到强迫对流的影响产生了自上而下的流动，这种流动导致了经由注水直接冲刷的正下方的格架内的水温极低，邻近套筒内的局部高温。

3) 入水管正下方的格架内的温度分布受到水流方向的影响，呈现出底部温度高、顶部温度低的分布特征。由于位于出口下方的套筒与邻近套筒组成了回路，使得周围邻近套筒的水温出现极大值。

4) 位于入水管正下方的贮存格架内的水流受到强迫对流的影响，是自上而下流动的。而远离入水管正下方的贮存格架内的水流受到自然对流的作用，均自下而上的流动。