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

(中国核动力研究设计院，四川 成都 610213)

反应堆控制棒驱动线作为反应堆内具有相对运动的设备单元。一般由燃料组件、导向组件、控制棒组件、控制棒驱动机构等部分组成。由于控制棒组件对堆芯反应性控制起决定性作用，因此一个合格的驱动线设计必须保证在紧急情况下，控制棒能够在规定时间内从反应堆内任意位置迅速插入堆芯，迫使反应堆停堆，防止事态恶化。因此，一直以来落棒时间也是控制棒驱动线设计首要考虑的指标，同时也是核电站安全分析的重要参考。由于物理样机造价昂贵，因此大批学者为准确计算落棒时间进行了大量研究工作[1-7]，这些研究工作将控制棒驱动线简化为一维水力模型，取得了大量的成果。这些成果作为控制棒驱动线设计的重要参考已经在核行业取得广泛应用，并形成了许多专用软件，例如CIGAL 软件(法国)、DROP 软件(美国)等。

但是，由于这些程序采用面向过程的编程思路进行编制，当驱动线结构发生变化时往往需要修改大量代码，扩展性较差。且由于采用一维水力模型，很多细节，包括流道形状变化均依靠引入实验修正系数来进行计算，因此对实验依赖性较高。随着CFD计算水平的提高和计算机技术的进步，一些学者选择采用CFD动网格方法对控制棒落棒行为进行仿真分析，例如，肖聪等基于 FLUENT 动网格技术对单根圆柱形控制棒建立了三维流体仿真模型来进行落棒行为仿真分析。同样，基于流体动力学(CFD)的动网格技术，在新的堆型上模拟了十字形控制棒组件在控制棒导向组件内的运动行为[8-9]。但是由于三维模型计算量过大，该方法只对驱动线的局部流域进行了仿真建模分析。到目前为止还没有针对完整驱动线的有限元计算分析。

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

1 流体仿真模型

1.1 网格模型

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

驱动线流场分为控制棒流域和驱动杆流域两个部分，忽略星形架等结构，并对流体域流道作了适当简化。以图1所示的控制棒流域为例，控制棒与导向管之间的流体域由Interface交界面分成了静止区域和运动区域两个部分，其中在运动区域，控制棒前端和后端分别定义了两个刚性面，其边界类型为Interior，主要作用是与Interface交界面将控制棒包围起来，所包围区域的网格可以随着控制棒进行刚性移动，而刚性面之外的运动区域网格则采用动态网格层变的方法进行网格更新，分割因子设置为 0.4，合并因子设为 0.2。这样就保证了控制棒能够在导向管中沿轴向自由移动而不会产生负网格。

图1 网格模型Fig.1 Mesh model

1.2 计算模型

本文采用FLUENT 14.5版本进行仿真计算。由于采用二维轴对称模型，因此计算时求解的是二维轴对称的连续性方程和动量方程，湍流模型采用k-epsilon湍流模型，使用SIMPLE 算法进行求解。

固体壁面均设置为固壁边界，其边界设置为无滑移条件，近壁面应用标准壁面函数。

压力插值格式选择标准格式，动量选择二阶迎风格式，湍动能及湍流耗散率皆选取一阶迎风格式。

反应堆冷态工况下，温度为 20 ℃，压力为101 kPa。整个流域与外界导通的部分，包括缓冲段的流水孔等，边界条件均设置为压力出口，表压为0。

2 分流域耦合方法

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

基本计算流程如图2所示，计算过程如下：

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

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

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

Δv=(M×g+F_d+F_c)×Δt/M

(1)

式中：M——控制棒和驱动杆的总重;

g——重力加速度。

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

v=v+Δv

(2)

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

5)驱动杆流域和控制棒单棒流体域利用DEFINE_CG_MOTION宏更新值更新Interior边界面(刚性面)的速度。

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

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

3 方法拓展性

本分流域耦合方法可以很容易由二维轴对称模型拓展到二维全流域模型甚至三维模型。甚至可以实现三维流域与二维流域之间的跨维度耦合。以三维模型为例，UDF代码部分几乎不需要进行修改。图1中的Interior刚性面在三维模型中替换为圆面，Interface交界面则为圆柱面。但是为了保证动网格动态网格层变方法顺利更新网格，需要对运动区域的网格沿轴线方向划分规整的分层网格。

4 仿真结果及分析

为保证计算结果的准确性，对模型进行了网格无关性研究，主要考虑控制棒间隙流域网格层数(r方向)以及流体域轴线方向(z方向)网格层数对计算的影响。如表1所示，为控制棒(低位)在导向管流速2 m/s条件下计算得到的稳态流体阻力值。从表中可以看出，控制棒流场计算对径向网格数和轴向网格数均比较敏感，随着径向网格层数和轴向网格层数的增加，计算误差在降低。综合考虑计算精度和稳定性，本文最终选取间隙流域网格层数(r方向)5层，轴线方向(z方向)网格层数6 000层作为网格划分方案。

表1 网格敏感性分析

一般而言，当驱动杆和控制棒从一定高度由静止状态在重力作用下沿着驱动线流道下落，初始阶段，由于流体阻力较小，落棒速度不断增加。此时驱动杆上端由于耐压壳内由于液体体积增加，造成了负压，并且导致流体沿着驱动杆和耐压壳之间的环形间隙向上流动，弥补驱动杆上端移动所形成的空腔。随着驱动杆速度的不断增加，驱动杆上端负压和环形间隙中的流体速度也不断增加，导致驱动杆流体阻力不断变大。这与图3(b)中0～0.7 s的变化趋势一致。

图3 流体阻力Fig. 3 Fluid resistance

对于控制棒而言，控制棒开始下落后，导向管中一部分流体通过导向管侧壁流水孔和底部端塞排水孔向外排出，另外一部分流体通过控制棒与导向管间的环形间隙向上流动，形成间隙射流；同时导向管内压强随着控制棒速度的不断增加而不断增大，这导致控制棒的流体阻力也在不断增加，但单根控制棒的流体阻力相对较小。当控制棒最终运动到缓冲段后，导向管侧壁流水孔被迅速遮蔽，排水作用减弱。控制棒与导向管间的环形间隙由于缓冲段缩口的存在瞬间减小。同时流通截面积突变和驱动棒对导向管底部流体的迅速挤压导致导向管内压强急剧增加，在控制棒底部形成很大的压差阻力，如图3(a)所示t=0.8 s左右形成的水力阻力脉冲所示。在此压差阻力的作用下，控制棒组件速度迅速减小。随着控制棒运动速度的减小，驱动杆和导向管内压差也相应不断减小，因此流体阻力也减小，并在低棒位附近最终趋于平缓。

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

图4 控制棒驱动线落棒速度和位移变化曲线Fig. 4 Variation curve of falling speed anddisplacement of the control rod drive line

图5 控制棒前端速度场Fig. 5 Velocity field at the front end of the control rod

图6 控制棒前端压力场Fig. 6 Pressure field at the front of the control rod

图5展示了控制棒下落过程中，控制棒前端速度场分布情况，图6则展示了控制棒前端压力场分布情况。可以看到随着控制棒的快速下插，控制棒前端形成较高的压力，导向管中的流体通过控制棒与导向管间的环形间隙向上流动，形成速度较大的间隙射流，这与前文描述的变化趋势相一致。

图7 压力分布曲线Fig. 7 Pressure distribution curve

图7给出了控制棒和驱动杆运动到某一位置时的表面压力，横坐标代表控制棒和驱动杆长度方向，其中0 m处表示棒和杆的上端，可以看到控制棒和驱动杆的上端均形成了负压，而在前段则形成正压，这与前文描述的规律一致。其中，驱动杆上端由于与耐压壳形成相对密闭的空间，因此负压作用明显，控制棒则由于处于开阔空间，因此负压作用较弱。

5 结论

本文提出了一种基于动网格的反应堆控制棒落棒行为分流域耦合仿真方法，该方法分别建立了控制棒单棒和驱动杆的二维轴对称模型，保证控制棒和驱动杆对应流域的网格能够根据落棒运动规律自适应地变化。两个流域在每个时间步长内交换流体阻力计算结果，并根据运动方程求解得到的速度来更新控制棒和驱动杆的运动状态，实现分流域耦合。该方法通用性好，计算中考虑了驱动线流道形状的影响，且在计算时间和求解精度之间取得了良好的折中，此外，该方法在计算条件允许的情况下，还能较容易地扩展到三维模型。本文采用该方法对某反应堆驱动线落棒行为进行了仿真，仿真结果表明，无论是速度、水力阻力、位移随时间变化，以及速度场、压力场分布情况等，均与一般理解的控制棒驱动线落棒规律一致，说明该方法计算结果可信。