一种基于重叠网格的反应堆控制棒落棒行为分流域耦合仿真方法

莫锦涛

（中国核动力研究设计院核反应堆系统设计技术重点实验室 成都 610213）

控制棒驱动线的落棒行为是指该设备在反应堆的寿命周期内，在高辐照、高温、高压并伴随有流致振动的条件下，控制棒组件（这里还包括驱动杆、可拆接头等部件）在驱动机构得到释放指令之后，在反应堆内的运动情况。由于控制棒驱动线落棒行为关系到事故工况下的安全停堆。因此落棒时间是反应堆结构设计时首要考虑的参数之一。

目前反应堆控制棒驱动线的落棒时间主要靠驱动线的相关综合试验来进行验证，但与此同时，控制棒驱动线的落棒时间的仿真研究也在开展。目前主要有两种方法用于反应堆控制棒驱动线的落棒时间的计算，一种方法是将控制棒驱动线简化为一维水力模型，对这种方法而言，需要求解的落棒动力学方程中无法通过理论计算得到水力阻力和机械摩擦阻力系数，这些数据只能依赖相关的试验获得，且无法考虑流道形状变化对仿真结果的影响。Donis等［1］首次通过求解一维动力学方程来计算得到控制棒落棒时间，在此基础上建立了数学模型用于分析单根控制棒落棒行为。美国、法国等为了对核电站压水堆控制棒落棒时间进行计算，均研制开发了相应的专用软件［2-4］，但由于这些软件均采用了一维动力学模型，引入了大量与试验相关的修正系数，因此在对不同堆型进行计算时需要对程序进行大量修改，拓展性较差。

另外一种方法是计算流体力学（Computational Fluid Dynamics，CFD）动网格计算方法。例如，肖聪等［5-6］基于计算流体力学动网格技术，对某反应堆的控制棒驱动线导向组件以及单根控制棒建立了相应的三维流体仿真模型，并利用FLUENT动网格技术进行瞬态动力学计算。但是由于三维模型计算量过大，因此该方法只对驱动线的局部流域进行了仿真建模分析，且由于计算量过大，计算时间过长，这样的方法距离实际工程应用还很遥远。到目前为止还没有针对完整驱动线的有限元计算分析。

随着CFD技术手段的进步，近年来，重叠网格技术取得了长足发展。动网格和滑移网格已经不再是模拟运动部件在流体中相对运动的唯一选择。重叠网格主要用来处理流体域中的部件运动，该技术的核心在于采用了背景网格和部件网格。部件网格用于模拟运动部件，通过部件网格和背景网格在不同相对位置处相互重合区域的网格进行数据交换，重叠网格较容易地解决了复杂相对运动下的流固耦合计算问题［7-15］。

本文提出了一种基于重叠网格的反应堆控制棒落棒行为分流域耦合仿真方法，该方法分别建立了控制棒单棒和驱动杆的二维轴对称重叠网格模型，保证控制棒和驱动杆的部件网格和背景网格能够根据落棒运动规律自适应地变化。两个流域在每个时间步长内交换流体阻力计算结果，并根据运动方程求解得到的速度来更新控制棒和驱动杆的运动状态，实现分流域耦合。本方法建模简单，方法通用性好，计算中考虑了驱动线流道形状的影响，且在计算时间和求解精度之间取得了良好的折中，此外，本方法在计算条件允许的情况下，还能较容易地扩展到三维模型。

1 流体仿真模型

1.1 网格模型

由于反应堆控制棒驱动线长度较长，而控制棒和驱动杆周围的流体间隙相对小很多，要获得与轴向长度和径向间隙长度相适应的网格尺寸，网格量将相当大。因此，为了减少网格模型大小，本文选用二维轴对称模型来对驱动线流场进行模拟，且为了保证计算求解精度和动网格的瞬时更新，整个流域采用四边形网格进行网格划分。

驱动线流场分为控制棒流域和驱动杆流域两个部分，忽略星形架等结构，并对流体域流道作了适当简化。以图1所示的控制棒流域为例，导向管内的流域被划分为背景网格，而围绕控制棒的部分流域被划分为部件网格，部件网格外边界定义为Overset边界。背景网格在对应控制棒部件网格径向间隙的部位进行了局部加密（本例中设置了10层网格），需要说明的是，控制棒部件网格径向间隙尺寸是由控制棒与导向管之间的最小间隙决定的（通常为缓冲段），而背景网格需要在控制棒部件网格径向间隙加密区扫过的区域进行相应的网格加密（如图1中10层网格深色区域所示），保证背景网格和部件网格的重叠区域至少有5层以上的网格。

图1 分流域耦合网格模型Fig.1 Multi-field coupling mesh model

1.2 计算模型

由于采用二维轴对称模型，因此计算时求解的是二维轴对称形式下的动量方程和连续性方程，湍流模型采用k-ε湍流模型。

本文计算条件选取反应堆冷态工况，常温常压，并假设控制棒从最高棒位由静止状态开始下落。计算流体域所有出口均设置为全压边界，流体域出口包括缓冲段流水孔、底部端塞排水孔等。固体壁面采用了无滑移无渗透的壁面条件并且采用壁面函数模拟边界层。模型求解器选用Coupled算法进行求解。为了兼顾计算稳定性和计算求解速度，除了压力和动量方程采用二阶离散格式，其他方程均选取一阶迎风格式。

2 分流域耦合方法

由于网格量和轴对称条件的限制，因此将控制棒驱动线的流场分为控制棒单棒和驱动杆两个部分。两个流场迭代计算同时进行，在每个时间步内进行数据交换以保持同步。

基本的耦合计算流程如图2所示，完整的计算过程如下：

1）控制棒单棒流体域计算为主计算流程，在每个时间步长内首先计算当前时间步长Δt，通过DEFINE_DELTAT宏修改时间步长，并输出Δt到文件time.txt，提供给驱动杆流域计算，驱动杆流域读入time.txt后同样通过DEFINE_DELTAT宏修改时间步长；

2）驱动杆流域积分计算驱动杆流体阻力，输出F_d到文件Force.txt文件，并等待更新速度文件；

3）控制棒单棒流体域读入文件Force.txt获得F_d，并积分计算控制棒流体阻力F_c（该值为单棒流体阻力诚意控制棒数），通过如下方程求解速度改变量：

同时更新速度，并写入文件vel.txt，根据更新速度值利用DEFINE_ZONE_MOTION宏更新动网格区域网格。

4）驱动杆流域读入vel.txt，得到速度值，利用DEFINE_ZONE_MOTION宏更新动网格区域网格。

5）数据交互完成并更新网格之后，两个流域同时开始迭代计算，迭代完成后进入下一时间步直到计算结束。

图2 分流域耦合方法流程图Fig.2 The flow chart of multi-field coupling method

3 仿真结果及分析

为保证计算结果的准确性，对模型进行了网格无关性研究，主要考虑控制棒部件网格和背景网格的重叠区域网格层数（r方向）以及流体域轴线方向（z方向）网格层数对计算的影响。如表1所示，为控制棒（低位）在导向管流速2 m·s-1条件下计算得到的稳态流体阻力值。表1中选取不同的径向（r）和轴向（z）网格数时的计算结果。从表1中可以看出，控制棒流场计算对径向网格数和轴向网格数均比较敏感，随着径向网格层数和轴向网格层数的增加，流体阻力值呈增大趋势，但计算误差的减小幅度在降低。综合考虑计算精度和稳定性，本文最终选取重叠区域网格层数（r方向）10层，轴线方向（z方向）网格层数6 000层作为网格划分方案。

一般而言，控制棒下落的过程中受到的力包括机械摩擦力、浮力、重力、流体阻力，在这些力的共同作用下完成落棒运动，为了简化，假设控制棒在对中条件下落棒，因此机械摩擦力可以忽略。控制棒在初始落棒位置时，由于落棒速度小，流体阻力很小，因此控制棒在重力作用下加速运动。在控制棒插入导向管过程中，不断对导向管中的流体进行挤压。受挤压的流体一部分被压出底部端塞排水孔和导向管侧壁流水孔，另外一部分则沿着控制棒与导向管之间的环形空间流出导向管。控制棒所受的流体阻力与控制棒端部的压差阻力以及控制棒与导向管之间的间隙大小有关，当间隙较小时，控制棒与导线管之间的速度梯度较大，流体作用在控制棒上的剪切力也较大。因此当控制棒下落到缓冲段时，带来流通截面积的瞬间变化，因此导向管内产生巨大压差。同时随着导向管与控制棒间的环腔体积瞬间减小，流体剪切力迅速增大，控制棒阻力也相应变大。在这些阻力的共同作用下，控制棒运动速度迅速减小，最终趋于平稳达到最终位置。

图3为控制棒驱动线落棒速度和位移变化曲线，结合图4中的控制棒和驱动杆流体阻力变化。可以看到与实验得到的控制棒落棒运动规律相一致，在进入缓冲段之前，控制棒在重力作用下不断加速，但是由于控制棒速度变快，相应的流体阻力变大，因此速度增加的趋势逐渐变缓。当控制棒进入缓冲段之后，由于控制棒底部形成的压差阻力以及缓冲段环形小间隙带来的摩擦阻力急剧增加，使得控制棒运动速度呈现断崖式减小的趋势，最终下降到一定速度后，稳定下移直到最终位置。从计算得到的控制棒落棒运动规律，可以看到与实际落棒运动特征非常符合，但是实验得到的落棒速度比计算值偏小，这可能是由于本文模型假设控制棒沿着轴线方向在对中条件下落，并未考虑实际落棒过程中由于错对中等因素带来的偏心、机械摩擦等，导致计算得到的落棒阻力偏小，因此落棒速度也就偏大。图5给出了控制棒进入缓冲段前后流水孔附近的压力分布变化，从图5中可以非常直观地看到，在进入缓冲段之前，流水孔附近的流体由于控制棒落棒运动的挤压而压力增加，部分流体由于高压而沿流水孔流出。而当控制棒进入缓冲段，流水孔附近的流体失去了运动挤压而压力迅速降低。

表1 网格敏感性分析Table 1 Mesh-density sensitivity analysis

图3 控制棒驱动线落棒速度和位移变化曲线Fig.3 The velocity and displacement curves of the control rod driving line

图4 控制棒驱动线落棒流体阻力变化曲线Fig.4 The fluid resistance curve of the control rod driving line

为了进一步验证耦合方法的正确，同时也对缓冲段流水孔的作用做一些研究。图6和图7分别给出了在流水孔直径减小一半、甚至不开孔的条件下，控制棒的落棒速度、流体阻力的对比。可以看到，在不开孔的条件下，由于缺少流水孔对导向管中受挤压流体的泄流作用，导致控制棒所受的压差阻力变大，控制棒下落速度整体变低（相比开孔条件下最高速度偏小0.7 m·s-1左右），进而使得控制棒落棒时间延长（相比开孔条件下到达缓冲段的时间变慢0.5～0.6 s）。此外，由于不开孔条件下，控制棒进入缓冲段的速度较小，因此缓冲段带来的阻力增加幅度也较小。

驱动杆所受流体阻力主要与落棒速度有关，驱动杆下降过程中由于耐压壳内液体体积增加，造成了负压，并且导致流体沿着驱动杆和耐压壳之间的环形间隙向上流动，弥补驱动杆上端移动所形成的空腔。间隙的流体流动同样给驱动杆带来了流体阻力。因此，在不开孔条件下，驱动杆流体阻力由于整体落棒速度的降低而偏小。

从图7对比还可以看到，落棒时间和流体阻力对孔径变化的影响敏感程度较小，但是对是否开孔的影响非常敏感。因此，流水孔的设置对控制棒落棒来说不可或缺，但是其开孔大小则可以进行一定的优化。

图5 缓冲段流水孔压力分布变化(a)控制棒进入缓冲段前，(b)控制棒进入缓冲段后Fig.5 The pressure distribution of the holes in buffer section of control rod drive line(a)Before the buffer,(b)After the buffer

图6 控制棒驱动线落棒速度变化对比Fig.6 The comparison of velocity curves of the control rod driving line

4 结语

本文提出了一种基于重叠网格的反应堆控制棒落棒行为分流域耦合仿真方法，该方法分别建立了控制棒单棒和驱动杆的二维轴对称重叠网格模型，保证控制棒和驱动杆的部件网格和背景网格能够根据落棒运动规律自适应地变化。两个流域在每个时间步长内交换流体阻力计算结果，并根据运动方程求解得到的速度来更新控制棒和驱动杆的运动状态，实现分流域耦合。本方法建模简单，方法通用性好，计算中考虑了驱动线流道形状的影响，且在计算时间和求解精度之间取得了良好的折中，此外，本方法在计算条件允许的情况下，还能较容易地扩展到三维模型。本文采用该方法模拟了某反应堆驱动线控制棒的落棒行为，仿真结果得到的速度、位移随时间变化与实验值符合良好，计算得到的速度场分布情况等与一般理解上的控制棒驱动线落棒规律特性一致，说明该方法计算结果可信。此外，还对不同缓冲段流水孔大小情况下控制棒的落棒速度、流体阻力变化情况进行了对比分析，发现落棒时间和流体阻力对孔径变化的影响敏感程度较小，但是对是否开孔的影响非常敏感。流水孔的设置对控制棒落棒来说不可或缺，但是其开孔大小则可以进行一定的优化。

图7 控制棒驱动线落棒流体阻力变化对比Fig.7 The comparison of fluid resistance curves of the control rod driving line