基于FLUENT 的UDS和UDF功能的中子扩散耦合计算研究与快堆应用分析

张雪贝， 王 驰， 陈红丽

(中国科学技术大学物理学院， 合肥 230027)

1 引 言

传统上，反应堆设计和反应堆安全分析一般都使用最佳估算程序. 随着计算机性能的提高和并行计算的发展[1-2]，反应堆的高置信度模拟在反应堆设计、方案优化和安全分析的研究中开始受到重视，只有考虑了反应堆模拟中的多物理反馈，才能实现高置信度模拟，而中子扩散和热工水力耦合计算是多物理耦合计算的重要组成部分[3-8]. 反应堆实际的运行就是中子学和热工水力相互反馈的一个过程. 反应性的温度系数(燃料温度系数和慢化剂温度系数)是反应堆正常运行时实现反应性控制的重要因素[9]. 要实现反应堆运行和瞬态工况的准确的计算，必须要实时考虑燃料温度、慢化剂温度和密度对局部中子通量和系统反应性的影响.

计算流体力学(Computational Fluid Dynamics, CFD)程序FLUENT能够通过对质量连续方程，动量方程，能量守恒方程的耦合求解实现反应堆堆芯和燃料组件的精细模拟. 因此将CFD模型与中子动力学模型耦合来进行反应堆安全分析的方法受到很大的关注[10]. FLUENT中的UDS(User Defined Scalar)能够对一类扩散方程运用Fluent内的求解器进行求解，其在多物理场耦合计算已得到广泛应用. 西安交通大学基于FLUENT的UDF和UDS功能开发TASNAM程序[11]，主要用于熔盐堆的中子扩散数值计算. 海军工程大学基于商用软件CFX的用户接口编程功能，增加三维时空中子动力学模型，并与CFD热工水力学进行耦合，对压水堆稳态下局部三维流动行为和三维物理特性进行了数值模拟[12].

本文基于FLUENT的UDF和UDS功能，对中子扩散方程进行定义，利用Gambit进行精细建模并利用FLUENT内的求解器在同一套网格下实现中子扩散与热工水力的耦合计算. 迭代计算流程如图1所示. 为验证耦合计算方法的可行性和数据传递的正确性，对5x5压水堆组件模型[13]进行建模和耦合计算，并将计算结果与其它耦合程序计算结果进行对比. 然后通过对M2LFR-1000[14]热组件的建模和耦合计算，得到组件稳态条件下的中子通量分布和温度分布，将热工水力参数(燃料最高温度，包壳外表面最高温度)和子通道程序(KMC-SUB)[15]的计算结果进行比较.

图1 耦合计算流程图Fig.1 The flow chart of coupling calculation

2 计算方法和数学模型

2.1 FLUENT的UDS模块

FLUENT的UDS模块能够利用有限容积法对一类方程进行定义和运用其内在的求解器进行离散和求解，形式如式(1)所示：

(1)

FLUENT中方程各项定义如表1所示.

瞬态中子扩散方程与UDS定义中各项对应关系如式(2)所示，式(2)方程左端第一项为非稳态项，左端第二项为扩散项，方程右端为源项.

对于稳态计算，略去方程中的非稳态项.

方程如式(3)所示.

(2)

(3)

式(3)中φg(r)代表中子通量，单位为cm-2s-1，D代表中子扩散系数，单位为cm，∑f为宏观裂变截面，单位为cm-1，∑g′→g代表从g′群到g群的宏观转移截面，单位为cm-1，∑r为宏观消失截面，单位为cm-1，ν为每次裂变释放的中子数 ，χg为中子的裂变谱.

有效增殖系数由式(4)计算：

(4)

表1 Fluent中的UDS对方程的定义

表2 5×5组件材料热物性

2.2 热工水力计算模型

冷却剂区域能量方程如式(5)所示：

(ρcoolantcPcoolantU·T)=βcoolantU·P+

(5)

燃料区域导热方程如式(6)所示：

(6)

(7)

γ为每次裂变释放的能量.

包壳和气隙导热方程如式(8)，式(9)所示：

(8)

(9)

2.3 材料温度截面库

通过DRAGON程序在已知材料核子密度条件下，同时考虑温度对材料密度的影响，通过能群归并，计算得到燃料，气隙，包壳，冷却剂在400, 500, 600, 700, 800,900, 1 000, 1 100, 1 200, 1 300, 1 400, 1 500, 1 600 K处的宏观截面，通过插值得到材料在400～1 600 K范围内的连续温度截面分布[16]，建立材料宏观截面和温度的函数关系，简写为式(10)，然后将宏观截面温度分布函数通过UDF添加到FLUENT的求解器中，通过读取每次迭代后每个网格的新的温度分布更新每个网格的材料宏观截面.

∑=f(T)

(10)

3 耦合计算方法验证

3.1 模型介绍

本文运用基于FLUENT的UDF和UDS功

能得到中子扩散和热工水力耦合计算方法对5×5压水堆组件模型进行计算，材料热物性由表2给出，模型结构和参数由图2和表3给出，冷却剂密度和比热随温度变化由式(11)，式(12)给出.

图2 5×5组件径向和轴向结构Fig. 2 Radial and axial geometric structure of the assembly

表3组件尺寸和材料参数

Tab.3Sizeandmaterialparametersoftheassembly

参数半径燃料棒半径4.1 mm包壳内径4.2 mm包壳外径4.8 mm燃料棒间距12.5 mm燃料高度3 m底部反射层高度0.2 m顶部反射层高度0.2 m燃料UOX(2%, 4% enrichment)冷却剂Water, 1 000 ppm boron气隙Helium, 0.1 MPa包壳Zircaloy-2功率12.5 MW能群数2能群分界能0.625 eV

3.2 建模与网格划分

使用Gambit对5×5压水堆燃料组件模型进行建模和网格划分，燃料棒径向网格划分如图3所示，组件径向网格划分如图4所示，中子扩散和热工水力计算采用同一套网格，其中轴向网格尺寸为0.1 m，模型网格总量3 100 000，利用Gambit的网格质量检验工具对所画网格进行检测，EquiSize Skew在0～0.4的网格数占99.32%，湍流计算采用k-epsilon模型.

(11)

(12)

3.3 计算边界条件

耦合计算边界条件如表4所示.

表4 耦合计算边界条件

3.4 计算结果及分析

为验证网格无关性，本文分别计算了930万网格和310万网格两套模型，得到310万网格条件下能满足计算精度要求，组件有效增殖系数参考值为1.171 09[13]，计算得到的组件有效增殖系数为1.171. 图5给出中子扩散和热工水力计算都达到收敛时2%和4%富集度燃料棒中心的轴向功率密度分布，蓝色点和黑色点为其它耦合程序结果，红色线和绿色线为本文的计算结果，4%富集度燃料棒中心最大功率密度计算值为4.25×108W/m3，2%富集度燃料棒中心最大功率密度计算值2.50×108W/m3，4%富集度燃料棒中心最大功率密度参考值为4.178×108W/m3，2%富集度燃料棒中心最大功率密度参考值为2.45×108W/m3.

由图5可以看出，计算偏差主要在组件模型的进口与出口处，功率密度参考值在组件进出口处有微小的上升，这主要是由于上下反射层的影响，使得部分中子被反射进入燃料区，进而引起功率密度的微小上升，本文计算时由于考虑方便运用Gambit建模和网格划分，忽略了上下反射层的影响，所以燃料棒轴向功率密度在进出口处未出现微小的上升. 图6为燃料棒外径处的温度分布图，图7为燃料包壳内径处的温度分布图，图8为燃料棒包壳外径处的温度分布图，图9为燃料组件进出口处的冷却剂和燃料棒的温度分布图. 图10为相邻燃料棒中心和对角燃料棒中心沿轴向的冷却剂温度分布. 图11给出z=0.0 m,y=0.0 m处沿x方向的温度分布(z=0为对称轴)，计算得到的4%富集度燃料中心最高温度为1 506.97 K，2%富集度燃料中心最高温度为1 066.42 K，4%富集度燃料中心最高温度参考值为1 502.22 K，2%富集度燃料中心最高温度参考值为1 047.60 K. 通过对压水堆5×5组件模型的耦合计算验证了利用Fluent 的UDS和UDF功能实现中子扩散和热工水力耦合计算的可行性和数据传递的正确性.

图3 燃料棒径向网格划分图

图4 组件径向网格划分图Fig.4 Radial mesh of the assembly

图5 燃料棒中心的轴向功率密度分布Fig.5 Axial power density distribution of the fuel rod center

Fig.6 Temperature distribution of the fuel rod outer diameter

Fig.7 Temperature distribution of the fuel rod internal diameter

图8 燃料包壳外径处温度分布

Fig.8 Temperature distribution of the fuel cladding external diameter

图9 组件进出口温度分布

Fig.9 Temperature distribution of the inlet and outlet

图10 相邻燃料棒中心处和对角燃料棒中心处冷却剂温度轴向分布Fig.10 Axial coolant temperature distribution of the contiguous fuel rod center and diagonal fuel rod center

图11 z=0.0 m, y=0.0m处沿x方向的温度分布Fig.11 Temperature distribution in the x axis direction (z=0.0 m,y=0.0 m)

4 耦合方法应用

通过对5×5压水堆组件模型进行耦合计算，验证了利用FLUENT的UDF和UDS功能实现中子扩散与热工水力耦合计算是可行的，现利用该耦合方法计算M2LFR-1000反应堆热组件模型.

4.1 模型描述

M2LFR-1000是作者所在课题组设计的模块化铅冷快堆，组件和燃料棒结构由图12和图13给出.堆芯燃料组件中的燃料棒采用正三角矩阵排布，构成的棒束为正六边形，被厚度为4 mm的组件盒包裹在其中，燃料棒的间距为14 mm，每个燃料组件包含有169根燃料棒. 燃料组件盒的内对边距为 185 mm，外对边距为193 mm，组件中心距为198 mm. 燃料芯块中心有直径为 1.9 mm 的中心孔，可以在燃料棒相同线功率密度条件下降低芯块中心温度，提高堆芯安全裕度. 燃料芯块外径为 8.6 mm， 采用添加少量 MA 或不添加 MA 的 MOX燃料.

燃料芯块和包壳之间有 0.15 mm 的间隙，其内填充有～0.5 MPa 的 He，填充气体压力高于一回路的运行压力，这样做一方面可以改善燃料芯块和包壳之间的间隙导热并提供惰性环境；另一方面可防止包壳因外压与自身蠕变坍塌造成与燃料芯块的接触；此外，还可以检验包壳是否破损. 燃料包壳为厚度0.55 mm的FMS T91，整个燃料棒的径向直径为10 mm，燃料棒活性区长度为1 000 mm.

选择综合性能比较好的 FMS T91 作为其堆芯结构和包壳的材料，FMS T91 的核素组成如表5所示[17]，冷却剂为Pb.

堆芯入口冷却剂温度设为673.15 K，堆芯出口冷却剂温度设为753.15 K，在包括设计基准事故(DBA)在内的各种工况下，燃料芯块最高温度应低于2 946.15 K[18]. 正常情况下，包壳最高温度应低于823.15 K[19]，并有足够的安全裕度.

图12 组件截面图(单位:mm)

图13 燃料棒截面图Fig.13 Fuel rod cross section

4.2 建模和网格划分

选取组件的1/6进行建模计算，热组件功率为3.6 MW，使用Gambit对M2LFR-1000组件模型进行建模和网格划分，其中轴向网格尺寸为0.05 m，模型网格总量为3 400 000，燃料棒径向网格如图14所示，组件径向网格由图15所示，利用Gambit的网格质量检验工具对所画网格进行检测，Equi Size Skew在0～0.4的网格数占96.55%，说明网格划分质量很好，湍流计算采用k-epsilon模型.

表5 FMS T91核子密度

图14 M2LFR-1000燃料棒径向网格划分图

Fig.14 Radial mesh of the M2LFR-1000 fuel rod

4.3 计算边界条件和物性参数[18]

计算边界条件由表6给出.

MOX热导率:

kMOX=(1.0/2.85×0.03+0.035+

T×(0.286-0.715×0.03))+(6400/T2.5)

e(-16.35/T)

(13)

MOX密度:

(14)

图15 M2LFR-1000组件径向网格划分图

表6M2LFR-1000组件耦合计算边界条件

Tab.6BoundaryconditionsofthecouplingcalculationoftheM2LFR-1000assembly

物理量边界边界条件类型边界值温度(T)组件入口定值673 K中子通量(ϕ)组件入口外推边界边界处梯度组件出口外推边界边界处梯度压力(P)组件出口定值101 KPa速度(U)组件入口定值(0,0,1.66) m/s

T91热导率

kT91=32.196-6.534e(-T/546.792)

(15)

T91比热

CPT91=404.336+T(0.881+T(-2.04×

10-3+T×2.6755×10-6))

(16)

Pb热导率

kPb=9.2+0.011T

(17)

Pb密度

densityPb=11441-1.2795T

(18)

Pb比热

CPPb=176.2-4.923×10-2T+1.544×

10-5×T×T-1.524×10-6/(T×T)

(19)

Pb粘度

viscosityPb=4.55×10-4e1069/T

(20)

4.4 计算结果分析

图16 和图17 为组件出口处的快中子和热中子通量分布(未归一化)，由于燃料区的裂变，快中子集中在燃料区域，热中子集中在冷却剂区域；并且由于燃料区域温度变化较大，造成此区域材料宏观截面参数发生明显变化，使得快中子和热中子在燃料区有明显变化；但是在冷却剂区域的快中子和热中子通量无明显变化，主要是由于堆芯冷却剂温度673～753 K左右，冷却剂Pb在此温度范围内材料宏观截面变化不大，所以对冷却剂区域中子通量分布无明显影响. 图18为组件外边界处的温度分布；图19为燃料芯块中心温度的轴向分布，由图19可以看出燃料芯块中心最高温度出现在z=0.6 m处；图20为z=0.6 m处，沿y轴方向上的组件温度分布，由图20可以看出由于中心孔的存在，使得燃料棒中心温度分布变得平坦，减低燃料芯块的最高温度，提高堆芯的安全裕度. 图21为轴向包壳外表面温度分布. 耦合计算得到的包壳外表面最高温度为774.6 K，燃料芯块最高温度为1 645.92 K，子通道程序KMC-SUB计算得到的包壳外表面最高温度为777.28 K，燃料芯块最高温度为1 647.17 K.

图16 出口快中子通量分布

图17 出口热中子通量分布

图18 外边界处温度分布

图19 轴向燃料中心温度分布Fig.19 The fuel pellet centerline temperature distribution along Z axis direction

图20 组件y方向温度分布(z=0.6 m, x=0.0 m)Fig.20 Temperature distribution along the y axis direction (z=0.6 m, x=0.0 m)

图21 轴向包壳外表面温度分布Fig.21 The cladding outer surface temperature distribution along z axis direction

5 结 论

本文基于FLUENT的UDF和UDS功能，对中子扩散方程进行定义，对反应堆组件模型进行精细建模，利用FLUENT内基于有限容积法的求解器对中子扩散方程进行迭代求解，同时进行热工水力的计算，实现了在同一个求解器和同一套网格下的中子扩散和热工水力的耦合计算.

通过对5×5压水堆燃料组件模型的建模与耦合计算，得到燃料棒不同边界处(燃料外径，气隙外径，包壳外径)处的温度分布，以及组件的进出口温度. 稳态条件下的组件有效增殖系数(1.171)与参考值(1.171 09)符合较好；将燃料棒中心轴向功率密度分布，相邻燃料棒和对角燃料棒中心的轴向冷却剂温度分布，x方向(z=0.0 m，y=0.0 m)温度分布与其他耦合程序计算结果进行对比，验证了耦合计算方法的可行性和数据传递的正确性.

将该耦合方法应用到M2LFR-1000中，对M2LFR-1000的热组件进行精细建模和耦合计算，得到组件内的中子通量分布和温度分布，并与子通道程序KMC-SUB计算结果进行对比，燃料芯块最高温度偏差为1.25 K，包壳外表面最高温度偏差为2.68 K，得到其热工水力学参数都在设计限值范围内.

下一步将考虑添加反射层的建模，然后进行瞬态的中子扩散与热工水力耦合计算.