堆芯燃料组件流阻特性模块式与整体式三维数值模拟方案比较

2020-05-30杨红义周志伟王予烨

原子能科学技术 2020年5期

杨红义，周志伟，林超，李淞，王予烨

(中国原子能科学研究院反应堆工程技术研究部，北京 102413)

快堆组件中相同功能组件按照其具体形状及所承担的具体功能任务分为不同类型的组件。每类组件对应的功率有所差异，故在反应堆运行时，不同类型组件对应入口流量不同。若燃料组件分配到的流量偏大，则堆芯出口处冷却剂钠的温度偏低，进而影响汽轮机参数，使得快堆经济性变差；若组件分配到的流量偏小，可能导致燃料组件包壳温度升高甚至超过极限温度，威胁快堆的安全。因此，须对燃料组件进行合理的设计和计算使其能匹配不同流量的冷却剂。而组件的阻力特性是影响堆芯流量分配至关重要的因素，所以建立合理的阻力特性计算方案成为当下亟待研究的问题。

堆芯包含大量组件，每盒组件又包含数十根燃料棒，因此要对整个堆芯进行三维流动换热计算模型的建模是极其不现实的。为方便示范快堆设计，中国原子能科学研究院堆芯及一回路项目组开发了适用于示范快堆热工水力设计及优化的程序，该程序需获得单个组件的流阻特性，以准确预估堆芯的流量分配是否合理。因此单个组件流阻特性的预测精度将很大程度上决定堆芯流量分配预测的准确性。

丁振鑫[1]对CEFR的Ⅰ-Ⅱ型小栅板联箱的流动阻力水力特性进行了实验。冯预恒等[2-4]则对 CEFR Ⅰ-Ⅱ型小栅板联箱及其节流件进行了稳态模拟计算，给出了节流件内的压力和速度分布，并与文献[2]中的试验结果进行对比，计算结果验证了三维数值模拟结果的可靠性，证明其完全可在核反应堆设计(至少在热工流体领域)中使用。王雨等[5]利用CFX对七根式小栅板联箱中快堆燃料组件管脚段均匀开口和非均匀开口两种模式下的速度场和压力场分布进行了数值模拟。李淞等[6]对169棒束的快堆燃料组件稠密棒束进行了三维数值模拟，根据模拟结果可更准确地预测棒束通道内的流动情况，可为今后稠密棒束组件水力学设计和子通道内流量测量试验提供参考。国外对池式钠冷快堆的流动特性也进行了大量的实验和数值模拟研究[7-8]。Padmakumar等[9]使用缩比模型对印度原型快堆(PFBR)进行实验，给出了反应堆内冷却剂流动分析、温度分布和流动振动特性。Péniguel等[10]对钠冷快堆冷却剂的三维流动特性和耦合传热进行了分析，但同样缺乏对局部压力分布及整体压降的模拟计算。燃料组件压降主要由3部分构成：管脚段压降、六角管棒束区压降和操作头段压降。上述相关研究都是对燃料组件各部分进行单独研究，目前尚未发现关于燃料组件全组件的流阻特性计算。戚展飞等[11]对带孔洞的复杂套管在不同冷却剂流动方向下的三维管段单相流动进行了数值模拟，研究表明，管道局部流动截面相对于流动方向的改变方式和改变顺序，是产生不同局部压降的重要原因。各组件部分之间过渡段也存在可观局部压降，分开计算并不能表征这部分的压降，同时分开计算也忽略了各部分之间流动的互相作用。为得到更为精确的单组件流阻特性，对燃料组件模拟件进行全组件三维数值模拟研究至关重要。

本文选用燃料区某两类组件作为研究对象，分别对其进行整体式和模块式三维数值模拟，获得其流阻特性并与试验结果进行对比，通过分析两种方法的精确度并结合计算资源耗费，综合评估两种分析方法的优势。然后深入研究速度和温度对两类组件流阻特性的影响，并分析这些因素对两种方案预测精度的影响。

1 模型介绍

本文以堆芯燃料组件为研究对象，选取其中两类燃料组件进行流阻特性研究。燃料组件结构主要包括操作头、六角管棒束段和管脚。燃料组件的棒束区与操作头结构一致，通过改变管脚区域的结构来改变整个组件的流阻特性。两类燃料组件的具体结构如图1所示，其中二类燃料组件在管脚(二类管脚)区域设置节流件来增加管脚区域的流动损失。

图1 燃料组件结构Fig.1 Structure diagram of fuel assembly

2 试验方法

试验模拟件与组件原型1∶1几何相似，且加工精度也与设计施工图纸要求保持一致。对于燃料组件，要求确保试验组件在水作为工质时的雷诺数与钠作为工质时的雷诺数相等，这意味着试验组件在水中及液态钠中应保持流速与运动黏度相同，进而充分模拟示范快堆实际工况中钠的流动状况与阻力特性，使得水力试验结果可应用于示范快堆的各项设计和验证。

沿燃料组件模拟件长度方向上冷却剂钠平均温度TNa约为420 ℃。由于水的运动黏度ν水与该温度下钠的运动黏度ν钠相等，可得到试验组件对应水温为93.2 ℃。由于两种工质的雷诺数和运动黏度均相等，可得到二者具有相

同的流速，进一步可得通过试验组件所要求的水工质体积流量Q水应等于组件平均温度下钠工质体积流量Q钠，即：

lNa=l水，νNa=ν水

VNa=V水

Q水=AV水=AVNa=QNa

同样，对于试验组件，进口钠温TNa为358 ℃时，可得到试验组件对应水温为84 ℃。

试验条件下采用恒定温度的去离子水作为试验工质，无法完全模拟实际的变温过程，因此需要在两个温度下分别进行试验，并进行温度折算。首先，在84 ℃与93.2 ℃进行组件无管脚试验，可在消除棒束区温差所导致测量误差的同时获得最贴近实际情况的棒束段温度；然后，在84 ℃下开展带管脚的全组件试验，获得相应的全组件压降；最后，84 ℃对应的管脚压降与93.2 ℃对应的棒束段压降构成组件的整体压降。

3 计算方法

3.1 网格模型

根据结构的特性，将燃料组件分成4个区域进行网格划分：管脚(两类)、棒束及操作头。不同类型燃料组件的棒束和操作头结构相同，只有管脚不同。对各区域网格以压降为标准，进行网格无关性验证。其中棒束区由于网格规模较大，棒束根据绕丝单个绕距为单位具有可复制性，因此本文取1个绕距进行网格无关性验证，结果如图2所示。

根据网格无关性验证结果，最终各区域具体网格数列于表1，具体网格划分如图3所示。

图2 网格无关性验证Fig.2 Grid independence verification

表1 各区域网格数Table 1 Grid type and number of each assembly

3.2 数学模型

燃料区组件管脚位于燃料棒下方，并无加热源，考虑到液体介质为钠液，可认为是动力黏度μ为常数的绝热不可压缩流体。因此，流动计算的控制方程可简化为以下两类方程[12]。

连续方程：

图3 各区域网格Fig.3 Grid diagram of each assembly

N-S方程：

目前，对各流体微团的紊乱运动进行平均处理加湍流模型求解方式仍是全三维湍流N-S方程求解的主要方法。在湍流模式中雷诺平均法在工程中的应用非常广泛，其基本思想是以雷诺假设为基础，将湍流的物理量用平均值和脉动值的和来表示，湍流场被视为平均运动场和脉动运动场的叠加[13]。本文选用目前工程中应用最为广泛的两方程模型作为湍流模型，其中两方程模型中RNGk-ε模型来源于严格的统计技术，它和标准k-ε模型很相似，但做了部分改进，使得RNGk-ε模型较标准k-ε模型在更广泛的流动中有更高的可信度和精度。本文选用RNGk-ε双方程湍流模型使方程组封闭。RNGk-ε模型具体方程[14]如下：

采用有限体积法对控制方程在空间上进行离散，选用SIMPLE算法对流场的离散方程进行求解。

3.3 边界条件选取

1) 整体式与模块式方案对比

为保证模拟条件和试验一致，本文选取全组件试验条件进行仿真计算，以去除数据处理过程产生的误差。该条件下介质为84 ℃的去离子水，此时水的密度为969.4 kg/m3，动力黏度为3.393×10-4Pa·s。选取进口流量为53.8 m3/h的工况进行计算和详细对比，进口为质量流量入口，出口采用压力出口边界条件。

2) 速度对压降的影响

为保证计算结果的可靠性，需用试验结果进行对比。本文选用整体式方案，并参照一类组件的全组件试验条件，选取介质为84 ℃的去离子水，进口为质量流量入口，出口采用压力出口边界条件，研究速度对压降的影响。速度与相应的体积流量与实验中部分值对应(表2)。

3) 温度对压降的影响

为保证计算结果的可靠性，需用试验结果进行对比。计算参照无管脚试验，计算值包含棒束及操作头区域，选用模块式方案进行计算。选取进口流量为53.8 m3/h的工况为基准，分别选取84 ℃和93.2 ℃进行计算和详细对比，进口为质量流量入口，出口采用压力出口边界条件。84 ℃去离子水的密度为969.4 kg/m3，动力黏度为3.393×10-4Pa·s；93.2 ℃去离子水的密度为963.4 kg/m3，动力黏度为3.03×10-4Pa·s。

表2 边界条件Table 2 Boundary condition

3.4 计算方案

二类燃料组件只有管脚结构不同，其余部件相同。因此计算区域可分为4个模块：一类管脚、二类管脚、棒束及操作头。本文对两种方案进行对比，并用试验结果进行验证。

1) 模块式方案

模块式方案对每个模块进行单独计算，获得每个区域的压降，通过各部分压力损失叠加获得燃料组件整体压降，如图4所示。

图4 模块式方案Fig.4 Modular calculation scheme

模块式计算的计算量相对较小、周期短、灵活性高，针对不同类型燃料组件，棒束和操作头的结构相同，具有可复制性，不用重复计算，只需对不同类型管脚进行计算即可。但模块式计算忽略了各模块间流动的相互作用的影响。

2) 整体式方案

整体式方案将4个模块组合成2组，将各区域通过interface进行连接，最终合成整体进行计算(图5)。

图5 整体式方案Fig.5 Integral calculation scheme

整体式计算考虑各模块之间流动的相互作用的影响，能更为精确地捕捉流动损失，获得的整体压降也更为精确。但整体式计算的计算量大、周期长、灵活性较低。

4 计算结果分析

4.1 模块式与整体式三维数值模拟方案对比

1) 流阻特性

对模块式和整体式两种方案的整体压降进行统计，并与试验结果进行对比，结果列于表3。从表3可看出，整体式计算的结果更为接近试验值。其中，一类组件整体式相对误差较模块式高4.9%，二类组件高1.5%。

对各子区域的压降进行对比，可获得不同计算方法影响最大的区域(表4)。可看出，绝对误差最大的区域为棒束区，相对误差影响最大的为操作头区域，而影响最小的为管脚区域。

表3 总压降对比Table 3 Comparison of total pressure drop

表4 各区域差值Table 4 Difference of pressure drop in each region

注：1) 绝对误差为各区域压降差值占组件整体压降试验值的百分比

2) 相对误差为相对各区域压降差值占各自区域压降整体式计算的百分比

由此可推测，整体式计算法对上游区域的流态较对下游的影响大，而下游对上游几乎无影响，所以位于最上游的管脚区域的整体式和模块式计算的压降几乎一致。而位于最下游的操作头区域虽然绝对误差不大，但相较其自身压降的相对误差接近30%。

对存在压降差值较大的区域进行分析，图6为距棒束区进口不同高度截面的轴向速度分布。可看出，模块式进口由于给定了流量进口，因此轴向速度分布均匀，而整体式由于上游管脚区域的影响，轴向速度存在高低速区。这种趋势直至截面距离进口高度为0.05 m时，两种计算模式的速度分布才逐渐趋于一致。因此在z<0.05 m区域两部分流动损失不一致，整体式由于速度分布不均引起的损失更大。

图7为操作头中心截面流线图，可看出整体式由于上游影响，入口流量存在偏离流动，在操作头中间部分存在旋涡。但由于操作头本身压降占总体压降的比例较小，因此该部分损失占总压降不到1%，可忽略不计。

2) 计算资源对比

使用计算资源为32核、256内存计算机，最多可计算2亿量级网格量，计算收敛时间约为120 h。

不同燃料组件只有管脚类型不同，在相同流量情况下，模块式计算法可将2个方案分解为4个模块，将4个模块的压降进行组合，即可得到整体压降。4个模块的网格量总和小于2亿，可同时计算，计算时间约为120 h。

整体式计算时每个计算模型网格大于1.7亿，2个计算模型总网格量为3.4亿，需分2次进行，计算时间约为240 h。

图6 棒束区z截面轴向速度分布Fig.6 z-section axial velocity distribution in bundle region

图7 操作头x截面流线图Fig.7 Streamline diagram in x-section of operating head

针对本文研究对象，模块式计算较整体式计算可节省一半时间。考虑到后续研究，堆芯共6种类型燃料组件，模块式计算所需时间为整体式计算的1/6。

4.2 速度对压降的影响

用模块式和整体式分别对一类组件和二类组件进行计算，并与试验结果进行对比，获得速度对压降的影响(表5)。可看出，随着速度的上升，压降上升；速度越大，整体式及模块式计算结果与试验值间的误差越大。

表5 速度对压降的影响Table 5 Effect of velocity to pressure drop

4.3 温度对压降的影响

由于全组件试验只对84 ℃去离子水进行测试，而无管脚试验对84 ℃和93.2 ℃分别进行测试，为保证计算结果的可靠近，选取无管脚试验所包含区域为研究对象，即计算模型只包括棒束区和操作头两部分。模块式和整体式计算所得的压降与试验结果的对比列于表6。从表6可看出，由于温度越高去离子水的密度越低，在进口体积流量一定的情况下，温度越高压降越小。整体式预测压降的相对误差约为3%，模块式预测压降的相对误差约为9%。温度升高，预测精度小幅提升。