核能的发展与应用对于我国综合国力的提升有重要意义。民用核能主要应用于核反应堆中，对于大多数压水堆，堆芯为棒束通道结构。因此研究棒束通道内的热工水力特性对于彻底了解核反应堆的运行机理至关重要。若反应堆在运行过程中发生了事故，非能动安全系统就会发挥作用以确保反应堆的安全，因此棒束通道内的自然循环流动特性是值得研究的。由于棒束通道结构的特殊性，有许多学者通过实验、数值模拟进行研究^[1-5]。实验研究可以得到更为准确、更令人信服的实验结果，数值模拟研究可以节省大量的人力、物力、财力。但对于想要较为精准地通过数值模拟来研究自然循环条件下棒束通道的热工水力特性，就需要将包括棒束通道的整个自然循环回路都进行模拟，这对于三维CFD数值模拟来说需要的时间更长，对计算机的要求更高。为解决这一困难，许多学者对三维CFD数值模拟软件与一维程序耦合进行了研究，李伟等^[6]证明了三维CFD数值模拟软件FLUENT与热工水力分析系统程序RELAP5耦合研究可以很好地分析简单的单相流或复杂的两相流问题；Palazzi等^[7-8]通过模拟管道中的摩擦阻力系数证明了STAR-CCM+与RELAP5-3D耦合的可行性；Gilman^[9]通过编写用户自定义程序建立了新沸腾模型，该模型能够准确预测CFD模拟中过冷沸腾流体的温度和热流密度。王宁波等^[10]对5×5棒束通道内的压降特性进行了数值研究，但对于单一的CFD数值模拟，仅能模拟强迫循环条件下的流动特性，无法模拟自然循环条件下的热工水力特性。孔松等^[11]用一维分析程序RELAP5对小型核动力装置进行了自然循环运行特性分析，但由于一维程序的局限性，很难对结构复杂的部件(例如反应堆堆芯)进行分析。

为探究反应堆堆芯内的固有安全性，本文研究自然循环条件下棒束通道内的流动特性，将研究棒束通道的三维CFD数值模拟与除去棒束通道外的自然循环回路其他部分的一维程序模拟耦合计算，以期用经济的方法计算得到更为准确的数值模拟结果。通过这种方法既可以将棒束通道通过三维CFD数值计算较为精细的模拟出来，又可以考虑自然循环回路对棒束通道的影响。自然循环回路管道多为普通圆管，采用一维程序进行模拟也可以得到很好的计算结果。本文采用三维CFD数值模拟软件STAR-CCM+与一维用户程序(User Code)耦合以实现自然循环条件下棒束通道的数值模拟。

1 三维与一维耦合方法

为确保三维与一维程序耦合数值模拟结果的准确性，首先要进行实验研究，将得到的自然循环实验数据用来验证数值模拟结果准确与否。实验装置为包含3×3棒束通道的自然循环回路，如图 1所示。本文所述三维与一维程序耦合模拟皆基于该试验台。具体实验方法见文献[2]。

为模拟自然循环条件下棒束通道的流动特性，本文采用基于有限体积法的CFD软件STAR-CCM+ 13.04对实验段3×3棒束通道进行三维数值模拟，并利用编程语言C++自行编写除去棒束通道外的自然循环回路其他管道，再将2个部分通过STAR-CCM+的用户程序自定义接口交互，以实现自然循环回路的完整数值模拟。

1.1 三维CFD数值模拟

棒束通道实验段结构如图 2所示。棒束通道实验段长度为800 mm，利用STAR-CCM+内的3D-CAD模型可直接进行建模。

1.1.1 网格独立性分析

本文采用多面体网格实现棒束通道主体结构的网格划分，采用棱柱层网格模型应用在所有的壁面边界。图 3给出棒束通道内生成的三维网格。

为了验证数值计算的准确性，首先需进行网格无关性验证。在网格独立的条件下对数值模型进行验证。以入口速度v=2 m/s的工况为例，对棒束通道进行了网格考核。以棒束通道出口温度T_o、棒束通道出口速度v_o及棒束通道进出口压降ΔP作为监测量。选取网格数分别为868 147、1 192 331、1 213 874和1 671 941的4套网格，当网格数从868 147变化到1 671 941时，棒束通道的出口温度T_o、出口速度v_o和进出口压降ΔP随网格数量的变化趋势如表 1所示。从表 1中可以看出，模型2中T_o的改变量为0.003 3%，v_o的改变量为0.15%，ΔP的改变量为1.03%，故认为网格数在1 192 331已达到独立。本文计算中选取网格数为1 192 331。

1.1.2 模型验证

为保证数值模拟的准确性，本文对竖直条件下的棒束通道进行模型验证。在图 1所示的实验装置中进行绝对压力为0.4 MPa，棒束通道加热功率为18 kW/m²的实验，将得到的实验结果与相同条件下不同湍流模型的计算结果进行比较。为找到最适合棒束通道的湍流模型，本文尝试了6种模型，分别为标准(Wilcox) k-ω模型、标准k-ε模型、可实现的k-ε模型、可实现的k-ε两层模型、雷诺应力湍流模型与SST (Menter) k-ω模型。数值模拟过程中将棒束通道入口设为速度入口边界，出口设为压力出口边界，棒束壁面设为定热流量，外壁面设为绝热，由此得到棒束通道进出口压降ΔP随雷诺数Re的变化曲线如图 4所示。从图 4中可以看出SST (Menter) k-ω模型可以很好地模拟棒束通道内的压降特性。故本文将选用SST (Menter) k-ω模型用于棒束通道的数值模拟计算。

1.2 一维User Code耦合方法

本文采用自行编写的一维用户程序(User Code)来模拟除棒束通道外的回路部分，包括预热器、冷却器与连接管路，如图 1所示。首先编写一维用户自定义程序，得到稳态条件下系统回路的自然循环流量，和该自然循环流量对应的速度，进而将其设为CFD数值模拟的入口速度，最后模拟得到自然循环条件下棒束通道内的流动特性。

因此，编写一维程序的目的就是找到实验回路达到稳定状态的自然循环流量值，即驱动压头与总阻力压头相等时的流量值。设冷却器与棒束通道中心高度差为L，棒束通道出口至冷却器进口冷却剂的密度为ρ_h，冷却器出口到预热器进口冷却剂的密度为ρ_c，采用换热中心假设，则驱动压头为ΔP_d=(ρ_c-ρ_h)gL。总阻力压头ΔP_z包括沿程阻力与局部阻力。沿程阻力即为各个管路的沿程阻力，具体计算方法为:

$ \Delta P_{f}=\lambda \frac{l}{d} \frac{q^{2}}{2 \rho A^{2}} $

(1)

式中：λ为沿程阻力系数；l为管路长度；d为当量直径；ρ为流过该管路的冷却剂密度；A为管路流通面积；q为质量流量。其中沿程阻力系数λ的计算方法^[12]为：

$ \lambda=\left\{\begin{array}{l} 64 / R e, R e<2320 \\ a+b R e^{-c}, 2320<R e<4000 \\ 0.3164 / R e^{0.25}, 4000<R e<1 \times 10^{5} \\ 0.0032+0.221 R e^{-0.237}, R e>1 \times 10^{5} \end{array}\right. $

(2)

式中：a=0.094k^0.225+0.53k; b=88k^0.44; c=1.62k^0.134; k=ε/d是相对粗糙度。

局部阻力包括棒束通道、预热器、冷却器处的局部阻力，以及各管路之间连接处，管路与棒束通道、预热器、冷却器连接处的局部阻力。局部阻力中除却棒束通道内的阻力外，其他部分相对较小，可以忽略不计。那么，总阻力压头表示为ΔP_z=ΔP_f+ΔP-ρgH，ΔP为棒束通道进出口压降，ρ为棒束通道的定性密度，H=800 mm为棒束通道长度。

因此，当实验回路中驱动压头与总阻力压头相等时，即ΔP_d=ΔP_z时即可输出自然循环回路质量流量q。具体程序框图如图 5所示。

值得注意的是，将STAR-CCM+模拟的迭代次数输入到一维程序中是为了编写实现STAR-CCM+迭代次数超过一定步数(例如100步)后再进行上述循环过程的语句，目的是让CFD达到相对收敛后再进行自然循环计算，从而避免计算结果的发散。

1.2.1 用户库

从上文可知，一维用户程序需要将STAR-CCM+中的计算结果作为已知量，并且需要将用户程序计算得到的自然循环速度传递给STAR-CCM+，使之作为棒束通道入口速度值。一维用户程序与STAR-CCM+之间的数据传递需要通过事先把用户自编程序编译成链接库，然后作为用户库加载到STAR-CCM+中的方式来实现。

1) 用户函数接口。

本文采用C++来编写用户函数，以实现得到自然循环回路流量的目的。通过库注册函数(uclib.cpp)指明用户函数类型为标量场函数“ScalarFieldFunction”，以此在STAR-CCM+中棒束通道入口处通过设置场函数的方法来设置入口速度值。注册方法可以通过调用ucfunc来实现，具体代码为:

ucfunc(main, “ScalarFieldFunction”, “IterationVelocity”);

其中，main为用户函数主函数名称，IterationVelocity为显示在STAR-CCM+下拉列表中的函数名。

另外，可通过调用ucarg注册用户函数所需的来自STAR-CCM+的参数，包括棒束通道进出口温度、压降，入口速度，迭代次数，具体代码为:

ucarg(main, “Cell”, "＄tinReport", sizeof(CoordReal)); /*棒束通道入口温度, K*/

ucarg(main, “Cell”, "＄toutReport", sizeof(CoordReal)); /*棒束通道出口温度, K*/

ucarg(main, “Cell”, "＄pjReport", sizeof(CoordReal)); /*棒束通道压降, Pa*/

ucarg(main, “Cell”, "＄vReport", sizeof(CoordReal)); /*棒束通道入口速度, m/s*/

ucarg(main, “Cell”, "＄Iteration", sizeof(CoordReal)); /*迭代次数*/

其中，Cell表示对于网格单元场的参数类型；Iteration为STAR-CCM+中原有的场函数，表示为迭代次数；tin、tout、pj、v为在STAR-CCM+中生成的报告，后缀加上Report可以同样实现场函数的功能；size为变量组分表中元素所需的储存(以字节为单位)，此大小可用于确保用户函数的精度与STAR-CCM+的精度相匹配，所有的场函数参数都应设置为CoordReal，即表示为双精度型浮点数据，另外尽管迭代次数Iteration为整型数据，也同样需要设定为双精度型浮点数据。值得注意的是，必须按照上述变量在用户函数中的所需顺序调用ucarg。本文中主函数的形参列表为：

void main(CoordReal *result, int size, CoordReal *tin, CoordReal *tout, CoordReal *pj, CoordReal *v, CoordReal *Iter)

其中，result为用户函数的返回值的组分表，对于本文来说即为自然循环速度值; size为result组分表中的元素数量。

至此可以得到库注册函数uclib.cpp代码如下：

#include "uclib.h"

void uclib()

{/*这里为上文所述的注册用户函数*/}

在与uclib.cpp相同的目录中创建文件uclib.h，以声明uclib.cpp中使用的函数。uclib.h文件中定义了在用户函数中使用的变量和函数类型，它对于所有代码都一样，具体内容参见STAR-CCM+用户指南^[13]。

2) 创建用户库。

至此，用户库源码已生成完毕，包括uclib.h、uclib.cpp、main.cpp, 由于在一维程序编写过程中需要查物性参数，因此还需调用物性库函数wasp97.h、wasp97.cpp。将上述文件同STAR-CCM+安装目录中的UserFunctions.lib编译链接成为动态链接库即可应用在STAR-CCM+中实现一维程序与三维CFD的耦合计算。

下面说明编译方法。本文在Microsoft Visual C++ 2013上进行编译，且Windows仅支持64位版本。打开VS2013 x64本机工具命令提示，将工作目录定位到当前工作目录，并使用以下命令将源程序uclib.cpp、main.cpp、wasp97.cpp编译到对象文件中：

cl /MD /D_WINDOWS /DDOUBLE_PRECISION -c ***.cpp

即可将源代码编译成二进制目标文件uclib.obj, main.obj, wasp97.obj。

下面说明链接方法。在VS2013 x64本机工具命令提示中输入如下命令：

link -dll /out:IterationVelocity.dll uclib.obj main.obj wasp97.obj "F:\\Program Files\\CD-adapco\\13.04.011-R8\\STAR-CCM+13.04.011-R8\\star\\lib\\win64\\intel16.3-r8\\lib\\UserFunctions.lib"

链接成功后可在当前工作目录中得到动态链接库IterationVelocity.dll。将该用户库在STAR-CCM+模拟中加载，可发现场函数列表中新增场函数User IterationVelocity，将其设置为入口速度幅值即可实现一维User Code与三维CFD的耦合计算，研究自然循环条件下棒束通道的流动特性。

2 结果分析

通过上述方法可以对自然循环条件下的棒束通道进行数值模拟。当系统压力为0.3 MPa时，在实验操作中，可以通过提高预热器加热功率或棒束通道实验段加热功率来增加实验回路的自然循环流量。对于数值模拟而言，为得到与实验相对应的多组工况，提高预热器加热功率转化为增加棒束通道实验段入口温度，提高棒束通道实验段加热功率可以通过直接在模拟中设置棒束热流量来实现。

若保持棒束通道加热功率不变，逐渐提高预热器的加热功率，可以得到系统回路自然循环流量随棒束通道入口温度变化关系。为方便与数值模拟结果对比，图 6展示了回路自然循环速度与棒束通道入口温度关系曲线，其中模拟速度值为棒束通道进出口速度平均值。从图中可以看出采用三维CFD软件STAR-CCM+和一维自定义用户程序耦合模拟可以很好地预测棒束通道内的自然循环速度值，可以认为STAR-CCM+与一维User Code耦合能够研究棒束通道内的自然循环流动特性。针对图 1所示实验过程中得到引压管2、3间的压降值ΔP₂，同样在模拟中得到相同部位的压降，并与实验值ΔP₂进行对比，如图 6所示，可以看出二者符合较好，最大误差在0.5%以内。

同样，若保持预热器加热功率不变，逐渐提高棒束通道的加热功率，可以得到系统回路自然循环速度随棒束通道热流量的变化关系，如图 7所示。对于压降ΔP₂也采用与图 6同样的处理方法，对比结果如图 7所示，发现二者仍然符合较好，最大误差在0.5%以内。可以看出无论是改变预热器功率或改变棒束通道加热功率，STAR-CCM+与一维用户程序耦合都能很好地预测棒束通道内的自然循环流动特性。

3 结论

1) 在STAR-CCM+中对3×3棒束通道使用SST (Menter) K-Omega湍流模型进行数值模拟，自行编写C++一维用户程序求解系统回路其余部分的流动特性，再通过STAR-CCM+的用户接口实现了二者的耦合。

2) 在改变预热器加热功率和改变棒束通道加热功率2种情况下将数值模拟计算得到的自然循环冷却剂速度值和棒束通道内压降结果与相应的实验值进行对比，发现二者符合较好，可以说明STAR-CCM+与一维用户程序耦合能够很好的预测棒束通道的自然循环流动特性。

本文对一维用户程序与三维STAR-CCM+之间的数据交互方法进行了说明，该耦合方法和实现流程为后续研究奠定了基础。