自然循环瓣阀开启行为的分流域耦合仿真方法

2022-03-22丰收王福鹏孙汝雷岳芷廷田瑞峰高璞珍

哈尔滨工程大学学报 2022年12期

丰收，王福鹏，孙汝雷，岳芷廷，田瑞峰，高璞珍

(1.哈尔滨工程大学核科学与技术学院，黑龙江哈尔滨 150001； 2.中国原子能科学研究院，北京 102413)

随着国家能源结构的调整，以核能、光伏等为代表的新能源发电产业受到越来越多的关注。“燕龙” DHR-400反应堆是中核集团设计的一种深水池式低温供热反应堆，反应堆采用了大量的非能动设计，其中自然循环瓣阀是安装在堆芯外层用于事故时形成“堆芯—水池”自然循环流动的关键设备。因此其性能的优劣直接决定了反应堆能否安全稳定的运行。

我国于 1980 年开始进行核能供暖技术的研究[1]，目前已经完成了多型核能供暖堆的设计工作，部分堆型进行了建造且完成了试运行工作。张亚军等[2]介绍了2×200 MW低温核供热产业化示范工程的概况、研究进展，总结了核供热堆的主要技术特点，并给出社会经济效益分析和应用前景展望。郝文涛等[3]对“NHR200-Ⅱ”型低温供热堆进行了介绍，阐述了核能供暖项目的安全性、先进性和可行性。柯国土等[4]从“燕龙”自身技术特点出发，开展方案设计优化、采用完全非能动余热导出技术、制定适用于“燕龙”的法规标准及实验验证策划。岳芷廷等[5]检验了DHR-200池式低温供热堆自然循环瓣阀在全厂断电叠加紧急停堆系统事故(SBO-ATWS)的安全性，验证了DHR-200池式堆的固有安全性。Jaakko Leppänen[6]以低碳政策的推行为切入点，总结并介绍了多国核能供暖反应堆的研发历史和技术特点。Gransell等[7]介绍了由芬兰和瑞典联合开发的SECURE反应堆的系统设计，对其非能动安全系统的作用机制进行了分析。Ville Tulkki 等[8]介绍了一种可以进行区域供热和工业供汽的小型模块化反应堆。Zhang Yixuan等[9]总结了低温池式供热堆的发展背景，介绍了核供热技术、安全性和经济性的相关研究。

计算流体动力学经过多年的发展，目前对定常问题的解决已经非常成熟，但是关于某些非定常问题的处理依旧面临着诸多挑战。重叠网格便是一种与传统动网格不同的网格处理方式，是由Steger等[10]提出，其基本思想是将计算网格分成了背景网格和若干个部件网格。多套网格之间相互独立，彼此的相对运动不对网格本身造成影响，从而允许计算模型可以有相对复杂的运动。在每一时间步内通过重叠区域的相互插值进行数据的交换，从而完成整个流场的计算。Robert等[11]为重叠网格的使用者提供了一个集成几何模型前处理、表面网格划分、体网格生成的用户操作平台，不仅提高了重叠网格的生成效率，而且能够极大的提高复杂物体的贴体网格的生成效率。Wang[12]提出了一种更为简单有效的自动挖洞方法，此方法在洞边界节点处采用 ADT 技术大大地提高了网格装配效率。在国内，庞宇飞等[13]提出了一种交点判别法进行挖洞的网格搜索，并且通过翼身组合体的例子进行了验证。经过前人的不懈努力，重叠网格技术取得了飞速的发展，目前已有利用重叠网格进行潜艇水下航行、汽车行驶过程流场分析、物体自由落体运动等多类具有大尺度位移的非定常运动问题中。

目前，有关自然循环瓣阀的热工水力学研究还鲜见报道，因此本文提出一种基于重叠网格理论的池式低温供热堆自然循环瓣阀开启行为分流域耦合的仿真方法。利用该方法分别建立了自然循环瓣阀和反应堆及堆水池的的三维重叠网格模型。本文采用该方法模拟了自然循环瓣阀的开启行为，获得了开启过程的速度场、质量流量等流场数据。

1 数学模型及数值方法

在建立的数值水池和堆芯中进行自然循环瓣阀旋转开启的模拟过程中，控制方程由连续性方程和动量方程(N-S)组成。

假定流体不可压缩，对事故初期自然循环瓣阀的开启过程进行数值模拟，需要满足以下控制方程：

连续性方程为：

(1)

动量方程为：

(2)

式中：ρ为流体密度，kg/m3；ui、uj(i，j=1，2，3)为X、Y、Z方向的速度矢量，m/s；gi为重力加速度分量，m/s2；μ为粘性系数，Pa·s。

重叠网格分为贡献单元和差值点，依靠二者的相对位置，可以求得贡献点与待差值点间的插值系数。根据这些插值系数，我们可以得到被插值点的任意变量值。

(3)

式中：ξi、ηi、ζi为插值系数；φi为贡献点上的变量值。

2 计算模型及数值计算

2.1 计算模型

由于DHR-400反应堆是双轴对称的结构，整个反应堆包含4条完全相同的一回路。为了在完整模拟事故初期自然循环瓣阀附近的流场变化的同时减小计算网格量的大小，本文的计算模型建立了1/3的反应堆结构原型，模型中包含了2条完整的一个回路，一个完整的自然循环瓣阀。堆芯及水池总高度3 m，堆芯半径1.5 m，水池半径5.0 m，堆芯及水池的圆心角120°，自然循环瓣阀闭合和完全开启时分别与竖直方向成±22.5°。计算模型如图1所示。

图1 反应堆数值计算模型Fig.1 Reactor numerical calculation model

自然循环瓣阀的结构如图2所示。反应堆正常运行时，射流管线的射流冲击会给自然循环瓣阀的阀瓣提供一个正压力；堆芯内冷却剂的高流速会使得自然循环瓣阀堆芯侧压强低于水池侧，形成一个向内的压力。在瓣阀自身重力及外力的耦合作用下，会使得自然循环瓣阀压在阀座上，成闭合状态。当反应堆一回路主泵失效时，射流管线关闭，主泵惰转不足以维持一回路冷却剂以额定流量流动，因此在自然循环瓣阀处有向外流动的趋势，施加给阀瓣一个向外的推力。当瓣阀内外的压差产生的推力和重锤的重力达到自然循环瓣阀的开启力矩时，自然循环瓣阀打开，形成“堆芯—水池”自然循环流动通路。为了使加快计算速度，减少计算资源的浪费，对自然循环瓣阀的模型结构进行了相应的简化:

1)删减瓣阀支撑肋结构。

自然循环瓣阀的支撑肋结构对整个计算过程中流场的影响甚微，删减其对整体计算结果的影响可以忽略。

2)简化自然循环瓣阀重锤结构。

自然循环瓣阀的重锤结构会对重叠网格区域的网格划分及数据插值造成较大的麻烦，因此对其进行简化。建模过程中略去重锤结构，重锤的作用通过UDF来实现。

图2 自然循环瓣阀结构示意Fig.2 Schematic diagram of the structure of the natural circulation flap valve

经过计算模型的简化，在不影响计算精度的前提下，既加快了数值模型的整体计算速度，还获得了贴近真实情况的流场信息。

为了使重叠网格部分拥有更好的对应关系，使得计算更加准确、收敛性更快，在自然循环瓣阀转动路径区域进行影响主体(body of influence，BOI)加密，加密区域切面如图3所示。

图3 转动路径加密区域切面Fig.3 Rotation path encrypts area slice plots

将计算模型导入后，建立流场计算域。由于自然循环瓣阀开启过程中会导致网格发生较大的畸变，故采用重叠网格的方式进行计算模型的网格划分。从反应堆运行安全的角度出发，考虑工程实际中自然循环瓣阀的加工误差、间隙装配及机械振动，遂将自然循环瓣阀模型的阀瓣半径设置成略小于自然循环瓣阀流道的半径(径向方向上有5 mm的间隙)。由于重叠网格方法不能处理固体域，所以采用“包覆-掏空”方法对自然循环瓣阀进行模型的处理，在自然循环瓣阀表面包覆一层5 mm厚的流体域，瓣阀的实体部分掏空处理，自然循环瓣阀流道的固体壁掏空处理。将自然循环瓣阀流体域定义为部件网格，其余部分定义为背景网格，且为了保证计算的精度及部件网格的瞬时更新，整个计算域采用六面体核心网格进行划分。瓣阀的网格如图4所示。

图4 自然循环瓣阀网格Fig.4 Natural circulation flap valve grid

在运动及网格的处理上，采用区域运动的方式，当自然循环瓣阀在内外压差的作用下运动时，自然循环瓣阀与表面流体域网格一起转动而不发生相对位移，通过部件网格和背景网格在不同相对位置处的重合区域的网格进行数据交换，改变流场参数在计算网格中的分布，从而插值获得自然循环瓣阀发生位移后的流场参数与网格对应关系。

整个计算域由部件重叠域及背景流场域组成，为了更精细地捕捉自然循环瓣阀附近流动的物理特性，自然循环瓣阀流道表面第1层网格节点的无因次长度y+控制在45左右，网格总数为810万左右。装配后计算模型切面网格如图5所示。

图5 装配后的计算模型切面网格Fig.5 Calculated model slice mesh after assembly

2.2 边界条件

计算采用隐式有限体积法离散动量方程，利用耦合求解器在时域中求解，计及重力影响并采用标准大气压作为参考压力进行初始化。对于此计算模型关注的重点为自然循环瓣阀开启过程中堆芯内外流场的变化，对边界层附近的流场变化关注度不高，且堆芯内冷却剂的流速很高，故使用高Re湍流模型k-epsilon进行计算，对流项采用二阶迎风模式以提高精度。

2.3 计算过程

在相关参数及边界条件设置完毕后，进行稳态工况的初始化和计算，为计算域提供一个初始的计算数值场。在稳态工况计算收敛后，改为瞬态计算，同时为瓣阀计算域设置压差敏感开启的UDF，使用Couple方法进行计算，瞬态时间步长采用0.005 s。计算过程如图6所示。

图6 重叠网格计算流程Fig.6 Overlay grid computing flowchart

2.4 网格敏感性分析

为了确保计算结果的准确性，对计算模型的网格进行了无关性验证。在保证自然循环瓣阀流道表面第1层网格节点的无因次长度y+控制在45不变的情况下，主要考虑自然循环瓣阀部件网格和背景网格的重叠区域网格数量对计算结果的影响。

如图7所示，随着网格数量的增加，自然循环瓣阀流出的质量流量呈现上升的趋势。网格数量从619万增加到810万时自然循环瓣阀的质量流量增加15%左右。网格数量从810万增加到1 050万时自然循环瓣阀的质量流量只增加4%。在计算精度允许的前提下，考虑计算所需时间及所耗费资源问题，顾选取网格量为810万的网格模型进行计算分析。

图7 网格敏感性分析Fig.7 Mesh sensitivity analysis

3 计算结果及分析

当反应堆正常运行时，进入反应堆堆芯的冷却剂被加热后通过主泵将冷却剂从一回路导出堆芯进行后续的冷却及循环。在事故发生时，主泵失效，一回路冷却剂的流动失去动力源头，导致堆芯内的冷却剂流量迅速降低。与此同时，关闭射流管线阀门。在瓣阀内外压力、瓣阀自身重力及自然循环瓣阀重锤重力的耦合作用下，自然循环瓣阀会打开，形成堆芯和水池之间的自然循环冷却回路。针对自然循环瓣阀开启的物理过程，以瓣阀标准安装高度为例进行分析。

3.1 自然循环瓣阀开启过程的速度场分析

在针对自然循环瓣阀开启过程的计算分析中，堆芯内的速度分布直接反应了堆芯内部冷却剂的流动状态，对自然循环瓣阀的结构设计及布置位置有着重要的影响。

在计算模型的对称面建立一个ZX平面，获得此截面上速度场的瞬态变化分布图。

为了模拟加工误差、间隙装配和机械振动导的自然循环瓣阀出现的漏流现象，在建模过程中自然循环瓣阀与流道径向上存在一个5 mm的间隙。从速度云图中可以看出，在自然循环瓣阀完全闭合时在瓣阀与流道之间一个漏流，漏流量为4 kg/h。此时漏流从瓣阀与流道之间的缝隙成射流状喷出，射流方向是垂直于自然循环瓣阀流道向堆水池。如图8是自然循环瓣阀闭合时流场的初始状态。

图8 自然循环瓣阀闭合时的速度场Fig.8 Velocity field when the natural circulation flap valve is closed

当自然循环瓣阀的开度为22.5°时(图9)，瓣阀与堆芯上部连接处流体速度较附近区域的流体速度大，且经过此处的流体有一个向左下方流动的状态，说明此时堆芯内的部分流体在压差的作用下向下方的自然循环瓣阀流动，判定堆芯上部存在大角度回转流或有旋涡存在。当流体在自然循环瓣阀出口的位置时，由于瓣阀下部流动阻力及堆芯水池压力梯度的存在，使得自然循环瓣阀出口下部的流体有一段微小的向上流动的趋势。由于自然循环瓣阀及瓣阀流道的阻挡作用，使得此位置瓣阀外侧出现一个绕流，绕流的外侧速度最大。当自然循环瓣阀完全打开时(开合角度45°)，如图10，瓣阀的流量达到最大。

图9 自然循环瓣阀旋转22.5°时的速度场Fig.9 Velocity field at a natural circulation flap valve rotation of 22.5°

图10 自然循环瓣阀旋转45°时的速度场Fig.10 Velocity field when the natural circulation flap valve is rotated 45°

3.2 自然循环瓣阀开启过程的质量流量变化

观察自然循环瓣阀出口的质量流量随时间变化的云图可以看出：当自然循环瓣阀完全闭合时(图11)，堆芯内的流体将自然循环瓣阀通道分成2部分，小部分流体泄漏进入水池内，大部分流体从自然循环瓣阀流道的一侧进入，另一侧流出。随着自然循环瓣阀的开启(图12～14)，从瓣阀流道流入堆芯的流体质量流量逐渐减小，流入堆水池的质量流量最大的位置逐渐移动到流道中心，流道中心位置流体的质量流量最大，周围的质量流量梯度变小。

图11 瓣阀完全闭合时流道入口的质量流量Fig.11 Mass flow of the runner inlet when the flap valve is fully closed

图12 t=0.12 s时瓣阀流道入口的质量流量Fig.12 Mass flow at t=0.12 s inlet of the flap valve runner

图13 t=0.2 s时瓣阀流道入口的质量流量Fig.13 Mass flow at t=0.2 s inlet of the flap valve runner

图14 瓣阀完全开启时的时流道入口的质量流量Fig.14 Mass flow at the runner inlet when the flap valve is fully opened

当事故发生时，主泵失效，自然循环瓣阀射流管线的流量为0。此时自然循环瓣阀在瓣阀内外压差、自然循环瓣阀重锤及瓣阀自身重力的耦合作用下获得了旋转角速度，瓣阀逐渐打开。随着瓣阀的开启，瓣阀流道的流量逐渐增加，从质量流量云图可以看出瓣阀流道内的质量流量逐渐增大，堆芯内的部分冷却剂从自然循环瓣阀处流出。

3.3 自然循环瓣阀的参数敏感性分析

为了探究自然循环瓣阀的结构参数对自然循环瓣阀开启行为造成的影响，本文针对自然循环瓣阀安装在400、600、1 000、1 400 mm的高度上进行了参数敏感性分析，获得了自然循环瓣阀在开启过程中的速度场、自然循环瓣阀及一回路的流量变化。图15所示为自然循环瓣阀的安装高度距离上升筒顶部400、600、1 000、1 400 mm时自然循环瓣阀流道的质量流量随时间的变化曲线。从曲线的变化趋势可以看出，不同安装高度的自然循环瓣阀在事故初期时从自然循环瓣阀漏流的流量基本相同，约为4 kg/h。随着自然循环瓣阀的开启，安装位置距离入口越近的自然循环瓣阀在同一时刻通过自然循环瓣阀流出的质量流量越多。同一时间不同安装高度的自然循环瓣阀质量流量如表1所示。当自然循环瓣阀的开启时间为0.075 s时，不同安装高度的自然循环瓣阀质量流量基本相同；当自然循环瓣阀的开启时间为0.3 s时，此时自然循环瓣阀处于全开的状态，全开的质量流量分别为311.0、320.1、364.6、387.2 kg/h。

图15 反应堆一回路质量流量随时间的变化曲线Fig.15 Curves of mass flow in the first circuit of the reactor over time

表1 不同工况的自然循环瓣阀在同一时刻的质量流量

如图16为不同的自然循环瓣阀安装高度工况下，自然循环瓣阀开启过程中堆芯一回路的质量流量变化曲线。从质量流量随时间的变化曲线得出：4种自然循环瓣阀的安装高度在事故时初期一回路的质量流量变化趋势及数值基本相同；不同安装高度的自然循环瓣阀的开启时间均为0.3 s左右；在自然循环瓣阀打开后的0.15～0.25 s内，自然循环瓣阀的安装高度对一回路的质量流量影响开始凸显，安装位置越高的自然循环瓣阀，其一回路的质量流量越大，结合自然循环瓣阀流道的质量流量变化可以看出，此时亦为自然循环瓣阀质量流量变化最明显的时间段。反应堆堆芯入口的质量流量随时间的变化曲线如图17所示。