钠冷快堆非能动停堆机构动导管共轭换热数值分析

任 逸，喻 宏

(中国原子能科学研究院，北京 102413)

钠冷快堆作为第4代核能系统中最成熟的快堆技术，在全世界已有超过420堆·年的经验[1]，我国已完成中国实验快堆(CEFR)的建造，目前正在进行中国示范快堆(CFR600)的设计。无保护失流事故(ULOF)是一种发生概率极小的严重事故工况，它是指失流事故时反应堆保护系统因机械故障或操作失误而不能正常工作，仅能依靠反应堆固有安全机制被动响应的严重事故工况[2]。为提高反应堆固有安全性，实现ULOF中控制反应性的基本安全功能，示范快堆设计了采用非能动机制的停堆系统，即液体悬浮式非能动停堆机构。

液体悬浮式非能动停堆机构由两部分构成，分别是非能动组件和非能动棒驱动机构。非能动组件由固定的组件套筒和可移动的非能动棒构成，在反应堆处于正常运行工况时，非能动棒悬浮于上工位；当发生失流事故时，非能动棒自动下落，插入堆芯实现停堆[3]。

动导管位于反应堆堆芯液态钠出口部位，其作用是为非能动棒驱动机构抓手提供运动通道和导向，液态钠自非能动组件出口部位沿动导管向上流动，最终流入冷钠池。

动导管内外两侧的环境差异较大，最大温差约为160 ℃。这样的温差所引起的热应力对动导管材料本身是一种考验，在长期运行时可能会使动导管寿命降低、断裂韧性下降。而动导管一旦断裂，将会出现严重的事故工况。这种在局部区域因较大温差导致的热疲劳问题已成为核电厂重点关注的安全问题。例如2008年10月，瑞典奥斯卡港核电厂(Oskarshamn Nuclear Power Plant)3号机组(沸水堆)停堆换料期间在控制棒焊接接头处发现了裂缝，其业主公司Tinoco等[4]进行广泛研究后，发现裂缝是由在控制棒和上部管道之间的环形间隙的冷层流(60 ℃)和旁路流动的热湍流(276 ℃)热混合导致的热疲劳产生的。

由于动导管两侧的热边界条件无法预先给定，而是受到液态钠和壁面之间相互作用的制约，因此动导管和内外两侧的液态钠构成一个典型的共轭换热问题，也称为耦合传热问题[5]。本工作采用开源CFD软件OpenFOAM中的chtMultiRegionSimpleFoam求解器求解该共轭换热问题。

chtMultiRegionSimpleFoam是OpenFOAM中较成熟且经验证的稳态共轭换热求解器，已成功应用于核能安全分析领域。2013年瑞典皇家理工学院(KTH)Ignacio[6]基于该求解器，对导致瑞典奥斯卡港核电厂3号机组控制棒破裂的低频温度波动现象进行了数值模拟研究；2017年，比利时冯卡门流体动力学研究所(VKI)Koloszar等与比利时核能研究中心(SCK-CEN)Keijers等[7]基于该求解器共同开发了用于MYRRHA热工水力行为数值分析平台——MyrrhaFoam；2016年，意大利米兰理工大学核工程系Pini等[8]在该求解器的基础上添加了内热源项，开发了新的求解器chtSourceMultiRegionFoam，研究有分布式内热源存在的熔盐回路的自然循环问题，且建造了实验台架对计算结果进行了对比验证。

本文采用Trelis前处理软件对钠冷快堆非能动停堆机构动导管和其内外两侧的液态钠建模并生成全域高质量六面体网格，基于开源CFD软件OpenFOAM中的chtMultiRegionSimpleFoam求解器，解决动导管和其内外两侧的液态钠构成的共轭换热问题，得到对动导管数值分析的温度分布结果，为力学设计提供重要的输入条件。

1 定义计算域与网格划分

定义计算域与网格划分是CFD计算的第1步，特别是对多个液体、固体域组成的复杂共轭换热问题，识别、选取、定义合理的计算域并对应划分高质量的网格尤为重要。本工作采用Csimoft公司开发的先进前处理软件Trelis进行几何建模与网格划分。

1.1 定义计算域

图1为通过Trelis建立的几何模型，除简化少量用于加工目的的倒角以外，基本保持原设计尺寸和形状不变。在图1的基础上，需进一步确定并生成相应的流体和固体计算域。动导管作为中间固体域(定义为Middle)，其内侧流体域(定义为Inner)是从非能动组件操作头进入的液态钠，外侧流体域(定义为Outer)是其周围1圈6个燃料组件操作头所形成的流体域。外侧流体域的选择可满足动导管共轭换热问题的研究和求解要求，同时可减少网格数、降低网格复杂度，提高求解效率。最终形成Inner-Middle-Outer 3部分构成的计算域。

图1 正常运行工况下动导管 及其周围组件操作头Fig.1 Guide tube and its surrounding assemblies’ operating heads under normal operation condition

Inner-Middle-Outer所构成的整体计算域的规模为：直径181.5 mm，高度580 mm。其中中间固体计算域Middle与图1所示动导管完全相同，对于流体计算域Inner和Outer，通过在对应圆柱区域中去除图1所示的固体部分得到，具体可通过Trelis中Boolean操作的Substract方法实现，生成的Inner和Outer计算域如图2所示。

图2 Inner(a)和Outer(b)流体域的生成Fig.2 Generation of fluid domains of Inner (a) and Outer (b)

1.2 网格划分

高质量网格的生成是CFD计算的前提条件，也是CFD工作中人工工作量最大的部分[9]，在本工作中约占总工作量的70%。本工作采用Trelis对整个计算域生成六面体网格。虽然六面体网格的生成耗时多、难度大，但其相比四面体网格有更好的精度、更高的效率。在六面体网格生成过程中，几何切割方案的制定是高质量网格生成的基础和难点，特别是OpenFOAM对网格质量的要求较高，如对网格的偏斜率、非正交度参数等指标有要求。

而要针对复杂几何沿流动方向进行“扫掠”生成高质量的全域结构化/非结构化六面体网格，必须对几何结构整体进行系统分析，从CFD计算的全局全过程考虑制定几何切割方案。例如图3所示Outer的切割方案中切割数约为300块。在此过程中Trelis和OpenFOAM需反复迭代，直到网格质量满足OpenFOAM的网格检查要求。

图3 Inner(a)、Middle(b)和Outer(c)切割方案Fig.3 Decomposition schemes of Inner (a), Middle (b) and Outer (c)

Inner-Middle-Outer生成六面体网格的数量和质量参数列于表1，其中网格已满足OpenFOAM内置的checkMesh网格质量检查工具的要求。由表1可看出，由于Outer区域几何形状复杂，切割难度大，使得其网格质量稍低于Inner和Middle。

2 计算模型

2.1 物理模型

一般的换热问题均会给定边界的温度分布、热流密度或这两者的线性组合，而共轭换热边界条件不能事先给定，如何构成定解问题并且合理地映射边界值是共轭换热问题的难点。chtMultiRegionSimpleFoam采用分区求解法求解共轭换热问题，即逐个求解各计算域的控制方程，再对相邻的边界的温度(T)和热流密度(q)耦合，重复此过程直到所有变量收敛[5]。

表1 六面体网格数量和质量Table 1 Quantity and quality of hexahedral mesh

流体域满足稳态不可压缩牛顿流体的控制方程。

连续性方程：

(1)

动量守恒方程：

(2)

能量守恒方程：

(3)

固体域满足传热方程：

(4)

式中：xi为方向坐标；ui为速度；ρ为流体密度；p为流体压强；μ为流体动力黏度；h为流体比焓；λl为流体导热系数；Tl为流体温度；Sh为流体内热源；λs为固体导热系数；Ts为固体温度。

在相邻计算域的交界面，对TW和qW进行耦合，例如对于Ⅰ区和Ⅱ区的边界需满足[5]：

TW|I=TW|II

(5)

qW|I=qW|II

(6)

此外，对于本工作中耦合边界面的处理有其特殊性，虽然Middle规则的几何形状可生成全域的结构化六面体网格，但由于Outer和Inner的不规则结构，如Inner在贯穿在斜面上的3个圆柱形流域，使得六面体网格生成时的扫掠方向受限，导致边界出现了不能保角映射的非结构化网格，如图4所示。针对此问题本文采用了OpenFOAM中的任意网格界面(AMI)技术。AMI技术可使每个面接受邻近的面和其交叉部分的权重的贡献，可确保在两个相邻边界面间插值且不会出现不连续问题。

图4 非保角映射边界网格示意图 (Inner-Middle-Outer)Fig.4 Boundary mesh diagram of nonconformal mapping (Inner-Middle-Outer)

2.2 湍流模型

雷诺时均(RANS)模拟方法中的k-ε两方程模型是工程中最常用的湍流模型之一，在液态金属管道流计算中已应用且经过实验验证。中国原子能科学研究院李淞等[10]应用标准k-ε模型对CEFR燃料组件管脚流量分配进行了数值模拟，且模拟结果得到实验验证。本文采用标准k-ε两方程模型[11]，湍流脉动动能k和湍流脉动动能的耗散率ε的稳态控制方程如下：

(7)

(8)

(9)

其中，u′i、u′j为速度脉动量。方程组包含5个可调整的常量，标准k-ε模型采用综合性数据得到的常数值可广泛用于各种湍流，Cμ=0.09，σk=1.00，σε=1.30，C1ε=1.44，C2ε=1.92。

此外，关于Inner和Outer的湍流边界条件，即湍流入口的k和ε，采用如下常用经验公式[11]计算：

(10)

I=0.16Re-1/8

(11)

(12)

l=0.07L

(13)

其中：Re为主流区的雷诺数；U为入口流速；I为湍流强度；L为当量直径。

2.3 物性参数

Inner和Outer计算域均为液态钠，其物性参数计算公式[12-13]参考国际原子能机构(IAEA)和比利时国家核能研究中心(SCK·CEN)发布的官方数据，具体如下所示。

液态钠密度：

ρ=[0.896 606 79+0.516 134 3×

(T×10-3)-1.829 721 8×(T×10-3)2+

2.201 624 7×(T×10-3)3-1.397 563 4×

(T×10-3)4+0.448 668 94×(T×10-3)5-

0.057 963 628×(T×10-3)6]×103

(14)

液态钠比定压热容：

cp=(38.12-0.69T-2-1.949 3×

10-2T+1.024×10-5T2)/22.99×1 000

(15)

液态钠动力黏度：

(16)

液态钠热导率：

λ=99.5-0.039 1T

(17)

Middle采用316不锈钢，其物性参数计算公式[12]如下。

316不锈钢密度：

ρ=8 084-0.420 9T-3.894×10-5T2

(18)

计算时认为动导管密度不变，取动导管平均温度477 ℃(750 K)时的值7 746 kg/m3作为动导管的密度。

316不锈钢比定压热容：

cp=462+0.134T

(19)

316不锈钢热导率：

λ=9.248+0.015 71T

(20)

2.4 边界条件

本工作计算的工况是正常运行工况，即非能动棒悬浮于上工位。处于正常运行工况的动导管内外组件操作头的结构和入口冷却剂参数分别如图1和图5所示。冷却剂入口边界值包括给定的温度值和流量值，设置为质量流量边界条件，该输入值取自设计值[14]，属于独立的边界条件；冷却剂出口设置为自由流出边界条件；固体部分均设置为固壁边界条件。

图5 正常运行工况下动导管 内外组件操作头入口冷却剂参数Fig.5 Coolant parameters for inlet of inner and outer assemblies’ operating heads of guide tube under normal operation condition

3 数值求解设置

3.1 离散格式设置

3.2 求解和算法设置

chtMultiRegionSimpleFoam求解器采用Simple算法求解流体域。因压力场在求解中收敛速度较慢[6]，故设置多重网格求解器求解压力场以加快收敛速度。为提高求解的稳定性，采用了亚松弛技术，松弛因子设置为0.3。此外，各场量内迭代求解的残差设置为1×10-7。

此外，本工作采用计算域分解法进行并行计算，OpenFOAM中的scotch分解方法是一种自动分解方法，其划分原则是使每块分割后的网格边界面最小化。本工作分别将各流体、固体域划分为6个计算，即6核并行计算。

4 计算结果和分析

4.1 网格无关性分析

网格无关性分析是CFD验证工作中最重要的活动之一[15]。为进行网格无关性分析，通过保持相同的几何切割方案和网格质量、仅改变网格划分间隔的方法形成由疏到密的3套网格，具体网格设置和最细网格尺寸列于表2。

Inner-Middle-Outer在3套网格下的比焓残差(相对残差)如图6所示。图中的横线为1×10-4，对于比焓残差，认为低于1×10-4即收敛。可看出，随着网格的加密，3套网格计算至5 000外迭代步均达到较好收敛。选取3套网格的5 000迭代步时Inner和Outer的出口流量及Middle的温度分布分别对流体域和固体域进行网格无关性分析。

表2 用于网格无关性计算的网格设置Table 2 Grid setting for grid independence calculation

图6 不同疏密网格下的Inner(a)、Middle(b)和Outer(c)比焓残差Fig.6 Specific enthalphy residuals of Inner (a), Middle (b) and Outer (c) under different grids

1) 流体域出口流量

不同疏密网格下Inner和Outer的出口流量列于表3。由表3可看出，随着网格的加密，Inner和Outer出口流量变化极小，可认为mesh3的流体域网格已满足网格无关性要求。

表3 不同网格下的出口流量比较Table 3 Comparison of outlet flows under different grids

2) 固体域温度场

Inner、Middle和Outer整体温度分布结果如图7所示。

由图7可看出，动导管外壁面的温差小于内壁面的温差，因此选取动导管温度场梯度大的内壁面，即Middle和Inner的交界面，高度选取4.38 cm，对不同疏密网格下围绕该高度一周的温度场计算结果进行比较，具体位置如图8a所示。由于Middle外部受到相对高温的Outer流体的加热、内部受到相对低温的Inner流体的冷却，特别是Inner中非能动组件操作头的3束分流对Middle的冷却，形成如图8b所示的波动的温度分布。可看出，从疏网格mesh1到密网格mesh3，其温度场分布差异趋于收敛，mesh3的固体域网格亦满足网格无关性要求。由图8可得到动导管在周向0.07 m范围内温差为97 ℃。

图7 Inner(a)、Middle(b)和Outer(c)温度分布Fig.7 Temperature distributions of Inner (a), Middle (b) and Outer (c)

图8 不同网格下动导管温度场比较Fig.8 Comparison of temperature distribution of guide tube under different grids

由以上对于不同疏密网格的流体域出口流量和固体域温度场的计算结果可得出，网格的划分方案满足网格无关性的要求，以下选择mesh3网格的计算结果进行分析。

4.2 内部流体域轴向温度分布

为验证选取的计算域高度合理，有必要对动导管内外两侧的流体轴向温度进行分析。本文选取动导管内部流体(Inner)轴向温度进行研究，结果如图9所示。由靠近出口处的h=420 mm和h=520 mm两图可看出，流体已经过较充分的混合，在出口的径向和轴向的温度变化均很小(温差约为20 ℃)，可判断选取的计算域高度合理。

4.3 动导管轴向温度分布

选取动导管3个轴向方向——Line A、Line B和Line C的温度分布如图10所示。由图10可看出，动导管在轴向0.17 m高度范围内的温差为80 ℃，且动导管轴向温度自下而上先急剧降低后再略有升高。前者是因动导管内非能动组件操纵头出口冷钠对动导管的冲刷冷却占主导作用，使动导管温度急剧降低；后者是由于动导管外部的燃料组件出口的热钠对动导管的传热占主导作用。

5 结语

本工作通过Trelis前处理软件生成全域高质量六面体网格，基于开源CFD软件OpenFOAM中的chtMultiRegionSimpleFoam求解器的应用，计算钠冷快堆非能动停堆棒驱动机构动导管和其内外两侧的液态钠构成的共轭换热问题，对网格无关性和计算域高度选取进行了验证和深入分析，得到了较为可靠的动导管数值分析的温度分布结果，得出动导管内壁面周向0.07 m范围内温差为97 ℃，动导管轴向0.17 m高度范围内的温差为80 ℃。在这样窄的范围内该温差将会产生较大的热应力，数值分析结果将为动导管力学设计提供重要的输入条件。

图9 动导管内部流体轴向温度分布Fig.9 Axial temperature distribution of fluid inside guide tube

图10 动导管轴向温度分布Fig.10 Axial temperature distribution of guide tube