基于CFD方法的控制棒下落行为研究

王振营，朱智强，张 蕊，刘志云，黄 宇

(1.中广核工程有限公司核电安全监控技术与装备国家重点实验室，广东 深圳 518172;2.武汉第二船舶设计研究所，武汉 430060;3.哈尔滨工程大学 核安全与仿真技术国防重点学科实验室，哈尔滨 150001)

控制棒组件是反应堆的控制部件，正常运行工况下具有实现反应堆启动、停堆和调节功率的功能.为了保证反应堆的安全性，要求控制棒驱动机构在正常工况下移动速度缓慢；在事故工况下，驱动机构能够自动脱开，控制棒组件在重力的作用下快速插入堆芯，整个过程可在2～3 s内完成.如果控制棒落棒速度过慢，或是紧急停堆时有一个反应性当量最高的组件对应的控制棒组件不能动作，将会导致非常严重的核事故[1-2]，因此控制棒下落时间是关系到核电厂安全的关键参数，对紧急情况下能否安全停堆至关重要，准确预测控制棒下落历程确定下落时间对反应堆的安全控制具有重要意义.

控制棒的下落时间与其所受到的水动力直接相关，不仅受其自身几何结构设计的影响，还受到冷却剂对它的作用力的影响[3-4].基于此，许多学者已开展了控制棒落棒相关的研究包括附加力场、常规作用力及结构优化等方面的问题[5-10]，但由于所处理问题的复杂性，往往采用实验方法拟合获得相应的计算关联式，然而实验关联式的应用范围受实验工况的限制，因此拓展性较差[11].

计算流体力学软件中的动网格技术被证明能够有效处理求解域及边界运动问题[12-13]，因此，本文旨在基于流体动力学软件FLUENT的动网格技术形成一套通用可靠的控制棒落棒数值分析方法，可用于多种堆型的控制棒落棒历程分析.

1 数学模型描述

1.1 基本守恒方程

本文采用的基本守恒方程如下：

连续性方程：

(1)

动量方程：

(2)

其中:ρ为液体的密度(kg/m3)，U为速度矢量(m/s).

基于湍动能k和耗散率ε的两方程模型是目前工程计算上应用最成熟、最广泛的湍流模型[14].本文选取的湍流模型为RNGk-ε模型，其表达形式如下：

湍动能方程：

(3)

湍流耗散率方程：

(4)

其中：μeff为有效黏度(Pa·s)；μt为湍流黏度(Pa·s)；p′为修正后的压强(Pa).

1.2 动网格基本方程

对于具有移动边界的任一控制体内通量φ的守恒方程通用表达式如下：

(5)

其中：ug为边界移动速度矢量，dV(t)为控制体V的边界，Γ是扩散系数(m2/s).

1.3 动网格技术

本研究中，由于控制棒及控制棒导向管的几何对称性，落棒过程中控制棒在水平方向的位移较小，在计算时可以忽略控制棒在横向的移动，仅考虑控制棒沿重力方向的位移，因此可以采用简化的一自由度(1DOF)求解器来求解该问题，并使用DEFINE_SDOF_PROPERTIES宏定义运动物体的质量和受力属性.需要注意的是，1DOF应用一般采用六面体网格并搭配使用纯铺层动网格模型.

2 计算模型和边界

2.1 几何模型

本文选取的实验是针对某小型反应堆的控制棒组件落棒历程设计的[15]，控制棒导向管为具有快插段和缓冲段的典型几何设计，快插段的尺寸比缓冲段大，控制棒在这段距离内不断加速，为防止燃料组件快速下落对燃料组件乃至堆内结构造成很大的冲击，在快插段的下部设置缓冲段来降低控制棒速度，缓冲段是在燃料组件导向管末端设置的内缩颈段，缓冲段上部及导向管下端通常布置有流水孔.实验对单个控制棒在导向管内的下落行为进行研究，控制棒尺寸与反应堆保持一致，单个控制棒的质量取4.2 kg，是将控制棒运动部件的质量平均分配到每根控制棒上获得的，图1为根据实验段建立的几何模型，具体的结构参数如表1所示.

图1 几何模型(横、纵向不同比例缩放)

表1 几何模型参数

2.2 网格划分

由于本文拟选择的控制棒落棒轨迹采用一自由度求解器，该求解策略需要采用六面体网格并搭配使用纯铺层动网格模型.而计算域属于典型的长纵横比几何，沿高度方向尺寸远大于径向方向几何尺寸，同时又含有局部复杂结构(流水孔和局部径缩结构)，对整体几何结构进行纯六面体网格划分难度极大且不可行，因此有必要采用混合网格，耦合使用结构化网格和非结构化网格，综合两种网格各自的优势，生成的网格如图2所示.

图2 网格划分

网格无关性检查结果如表2所示，分析不同网格数下流场参数，确定最理想的网格数为1.66×106.

表2 网格无关性检查

2.3 计算边界及求解设置

控制棒在重力的作用下不断向下运动，控制棒底部的网格不断被压缩消失，网格数目减小；控制棒顶部的铺层网格不断被拉伸分割为新的网格，网格数目增加，图3给出了落棒过程中动流体域内的网格变形情况.

图3 控制棒沿导向管下落网格变形情况

由于控制棒的下落特性不是预先给定的，需要通过受力分析确定，控制棒的位移主要体现在竖直方向，因此采取“1-DOF”模型来求解被动型动网格问题.实验在常压下进行，压力为0.1 MPa，温度为20 ℃.计算域出口边界条件设置为压力出口，压力为1.01×105Pa.控制棒组件初速度设置为 0 m/s.

3 计算结果及分析

图4给出了本文数值模拟获得的控制棒的速度-时间曲线.从图4中可以看出，控制棒在下落初期速度不断增加，当达到最大速度2.0 m/s左右时速度逐渐衰减，这是由于控制棒的质量是给定的，在下落的初始阶段，重力远大于控制棒所受的浮力和摩擦力，加速度的方向向下，控制棒不断加速.随着控制棒的速度逐渐增大，控制棒所受的浮力和摩擦力逐渐增大，加速度不断减小，控制棒的速度进一步增大.当控制棒的速度增大到最大值时，控制棒的加速度为0，重力与浮力、摩擦力共同达到平衡.控制棒进一步向下运动时，浮力、摩擦力等方向向上的合力大于重力，控制棒速度逐渐减小.控制棒进入缓冲段之后，由于缓冲段的直径小于快插段，控制棒所受到的摩擦力明显增大，竖直向上所受合力进一步增大，加速度增大，控制棒的速度迅速降低.

图4 控制棒速度-时间曲线

图5给出了计算得到的控制棒位移-时间曲线与实验测量之间的对比情况，二者符合良好.

图5 控制棒位移-时间曲线

图6所示为控制棒导向管内压力随时间变化云图.可以看出，随着控制棒不断进入导向管，导向管内压强不断增大.当控制棒继续向下运动进入缓冲段后，控制棒与导向管间的环形间隙的尺寸瞬间减小，形成很大的摩擦阻力.随着摩擦阻力的增大，控制棒的下落速度逐渐减小，导向管压强也逐渐减小.

图6 导向管压力云图

控制棒下落时间一直是关系反应堆安全的关键参数，通过数值模拟可以进一步将下落时间细分为快插时间、缓冲时间，分别是515 ms和140 ms，控制棒最大速度为2.16 m/s，表3将实验和数值模拟结果中这些关键参数进行对比，考虑到数值模拟几何与实验段尺寸之间的差异，两者之间的误差较小.

表3 数值模拟和实验数据对比

4 结 语

本研究基于计算流体力学软件中的动网格技术对控制棒下落行为进行数值模拟研究.计算获得单根控制棒落棒的位移、速度和时间以及控制棒导向管内的压强分布等结果，并基于公开发表的单根控制棒落棒实验数据对数值模拟结果进行验证，结果表明本文提出的数值模拟方法可以准确地预测控制棒下落过程中的位移及时间相关参数，能够为反应堆控制棒落棒相关结构优化设计和安全分析提供依据.