作为堆芯的入口边界，压力容器下腔室出口的冷却剂流量、温度及硼的分布对反应堆的安全性和经济性有着重要影响。因此，对压力容器中冷却剂流动特性进行研究成为反应堆安全分析的重要组成部分。随着CFD技术的快速发展，利用CFD分析方法对压力容器进行数值模拟，逐渐成为与实验同样重要的研究手段。中国核动力研究院的眭曦等^[1]基于下腔室实验对CFD计算中不同交混模型进行对比，认为组分输运模型可以较好的描述由于浓度梯度导致的硼扩散现象。郭超等^[2]基于CFX程序进行下腔室全尺寸计算，并对下腔室流场及压力场进行研究，验证下腔室流量分配设计合理性。Cho等^[3]基于ROCOM实验对CFD程序进行验证，验证结果表明CFD程序在适当的网格处理情况下可以合理的预测压力容器内的硼稀释现象。

大部分关于压力容器的实验的装置按比例对真实压力容器进行缩放，但其复杂的下腔室结构和依然较大的模拟域导致CFD模拟所需的网格数量庞大，对计算使用的硬件提出了较高的要求。同时，冷却剂在压力容器下降段与下腔室结构差异很大(下降段为环段结构，下腔室为球形且内含孔板结构)，使用同一湍流模型可能无法保证精度与效率的最佳平衡。为解决这些问题，本文将分段技术和MRANS方案应用在压力容器CFD计算中。Navarro等^[4]采用分段模拟技术对Karoutas等^[5]的5×5棒束实验进行模拟，以减少每次计算的计算负荷；Li等^[6]对分段模拟技术在棒束通道仿真中的应用特性进行了研究；Chen等^[7-8]针对带交混翼棒束通道不同区域流体流动状态不同的特点，基于分段模拟技术提出了多湍流模型(Multiple RANS models, MRANS)方案，在保证计算精度的同时提高计算效率，并提出了DP-CFD方案以解决堆芯大流域CFD仿真问题。值得注意的是，绝大部分的分段CFD分析技术和MRANS方案的讨论都是针对棒束几何模型，在其他几何模型中的适用性仍然需要讨论。

ROCOM实验^[9-10]被广泛应用在CFD计算方法的验证^[11]中，许多学者基于ROCOM实验装置进行了CFD分析，其中包括网格敏感性、湍流模型和边界条件效应等等^[12-14]。本文以ROCOM实验为对象，针对分段技术与MRANS方案展开研究，提出合理的分段计算与MRANS方案，减小计算负荷，提高计算效率，并保证CFD计算的准确性。

1 仿真对象介绍

本文以欧洲FLOMIX-R基准实验中^[9-10]的ROCOM实验装置为仿真对象，对稳态运行的压力容器的冷却剂交混现象进行研究。ROCOM实验装置以KONVOI三代压水堆为原型，按1∶ 5的比例建造，反应堆压力容器(RPV)为其主体设备，包括4个冷却剂入口、出口、下降段(downcomer)及下腔室(lowerplenum)，下腔室布置有空心孔眼滚筒，堆芯入口布置支撑板。装置具有4个完整的回路，并配有泵、阀门及蒸发器等设备，通过对泵的流量控制可以实现不同工况的实验。在进行混合实验时，实验装置在室温和环境压力下运行，示踪剂氯化钠溶液在入口上游注入主冷却剂系统，并通过装置进行混合，确保示踪剂均匀分布在压力容器入口横截面处，压力容器上布置有4处测点(容器入口、下降段上部、下降段下部及堆芯入口)对示踪剂进行跟踪，具体参数如图 1所示。

2 分段CFD技术研究 2.1 数值建模 2.1.1 几何模型和网格划分

根据Boumaza等^[14]的研究，环路对发生在RPV的交混现象的影响很小。因此，在不影响计算准确性的前提下减少计算所需网格，本文使用300 mm的圆管代替环路作为RPV入口。堆芯采用高度为800 mm的圆柱代替，圆柱以堆芯入口支撑板上直径30 mm的冷却剂流量分配孔为基准进行延长。

图 2展示了本文所使用的网格模型。如图 2所示，本文采用多面体网格对RPV模型入口环腔和下腔室进行网格划分，这是因为多面体网格与四面体网格相比，在达到相同精度的前提下具有更好的收敛性^[6]。多面体网格的生成方式采用将四面体在FLUENT软件中转化的方式生成。下降段则采用六面体网格。表 1展示了在网格敏感性分析中使用的不同划分方案的网格量。图 3显示了不同网格方案中堆芯入口处标量最大和平均值随时间的变化趋势。由于网格划分占据了前处理中的绝大部分时间，同时无止境的增加网格量会增加大量的计算时间。因此，尽管未达到所谓的网格无关解，考虑精度与计算时间的平衡的情况下，本文选择了方案2作为本文的计算网格划分方案。

2.1.2 边界条件和湍流模型

本文以ROCOM_STAT_02实验工况为仿真对象，其是针对低流速下的稳态混合进行的实验，实验过程中没有发生热传导且操作环境为室温和环境压力，因此选择25 ℃，1.01×10⁵ Pa的水作为计算工质，采用0.73 m/s(25%正常流速，即各环路的质量流率为46.25 m³/h)的均匀速度分布作为入口边界条件，出口采用压力出口边界条件。作为硼替代物的混合标量采用均匀分布的方式由入口1中注入，注入时间共35 s。本文采用的时间步长与实验中测点的时间监测间隔相同，均为0.05 s(频率20 Hz)。本文的分段计算方案研究中, 均使用RSM湍流模型。

2.2 分段方案研究

RPV根据流域特点可以分为：下降段(DC)和下腔室(LP)。在建模和网格划分过程中，将DC和LP从整体RPV中拆分，得到代表各区域的独立模型。本文利用FLUENT软件带有的界面数据读写功能，通过编写的Python控制程序控制DC和LP之间的数据传输，完成分段计算的自动计算流程。图 4展示了DC和LP的数据传递过程。

在文献[4, 6-7, 15]的分段计算中，将计算区域沿流动方向分为若干段，并将前一段的出口的速度场结果作为下一段的入口边界条件，如图 4中黑线的传输方案。这种计算方式相当于把流场下游对上游的影响切断，意味着假设了下游的几何与流动状态对上游没有影响。这种数据传输方式可以保证在没有回流模型中速度场的准确传播，但由于缺乏下游模型，上游的出口压力场将无法准确得到。综上所述，由于缺少数据通信，会导致上游出口面与下游入口面的压力分布不一致，无疑会引入一定的误差。

针对上、下游交界面压力不一致的问题，本文提出了一种基于区域重叠的双向数据传递分段计算方案。区域重叠方案如图 5所示，DC域向下拉伸了150 mm的圆环柱与LP域重叠。在计算过程中，DC将LP的入口所在面的速度和其他标量数据作为入口边界条件传递至LP，而LP将DC的出口所在面的压力分布数据作为出口边界条件传递至DC，其与文献[4, 6-7]的单向数据传输的不同之处在于数据传输包括了图 4中红线的传输内容(压力场)。图 6展示了不同的分段方案交界面处的压力对比图。

如图 6所示，相较于pressure-outlet出口边界的均匀压力场，outflow出口边界的压力分布在趋势上与整体计算相似，这与Li等^[6]所得到结果相同。然而，由于单向数据传输的自身缺陷，使用outflow边界在交界面处的压差与整体计算依然有一定差距(整体计算中界面压力最高点与压力最低点的压差为100 Pa，而outflow计算为40 Pa)。而双向数据传递的分段方案考虑了上下游的相互反馈，从图 6可以看出，双向数据传递的分段方案与整体计算最接近。

图 7、8分别为不同的分段方案与整体计算的结果对比。图 7为DC下降段下部的3个测点的速度。考虑到计算几何具有对称性, 3个测点选取在同一个象限内。图 8为LP中16号通道与42号通道堆芯入口平面的速度。选取这两个通道是因为两个通道分别在下腔室空心孔眼滚筒的内外两侧，具有一定的代表性。从图中可以看出，采用双向数据传递的分段方案与整体计算结果对比具有较好的一致性。

3 多湍流模型方案研究 3.1 MRANS方案设计

受到计算资源的限制，DNS和LES模型在工业CFD仿真领域使用较少，通常采用基于雷诺时均方法的RANS模型进行计算。雷诺时均的N-S(Navier-Stokes)方程可以表示如下：

$ \frac{\partial \rho}{\partial t}+\frac{\partial}{\partial x_{i}}\left(\rho u_{i}\right)=0 $

(1)

$ \frac{\partial}{\partial t}\left(\rho u_{i}\right)+\frac{\partial}{\partial x_{j}}\left(\rho u_{i} u_{j}\right)=-\frac{\partial p}{\partial x_{i}}+\tau_{i j}+\frac{\partial}{\partial x_{j}}\left(-\rho \overline{u_{i}^{\prime} u_{j}^{\prime}}\right) $

(2)

式中：ρ为流体密度, u_i为在i方向上的平均速度分量, p为压力, τ_ij为黏性应力张量。雷诺应力$-\rho \overline{u_{i}^{\prime} u_{j}^{\prime}}$的求解是RANS模型的关键。

标准k-ε模型(SKE)、realizable k-ε模型(RKE)、renormalization group k-ε模型(RNG)及shear stress transport k-ε模型(SST)等2方程模型基于Boussinesq假设和湍流黏度的概念，通过下式求解雷诺应力。

$ -\rho \overline{u_{i}^{\prime} u_{j}^{\prime}}=\mu_{t}\left(\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}\right)-\frac{2}{3}\left(\rho k+\mu_{t} \frac{\partial u_{k}}{\partial x_{k}}\right) \delta_{i j} $

(3)

式中：k为湍流动能, μ_t为湍流黏度。湍流黏度是湍流动能和湍流耗散率或比耗散率的函数，因此求解两方程模型的雷诺应力需要建立湍流动能输运方程和湍流耗散率或比耗散率输运方程。

雷诺应力模型(RSM)摒弃了各向同性的湍流黏度概念，直接构造了关于雷诺应力的输运模型，一般如下式所示。

$ \frac{\partial}{\partial t}\left(\rho \overline{u_{i}^{\prime} u_{j}^{\prime}}\right)+\frac{\partial}{\partial x_{k}}\left(\rho U_{k} \overline{u_{i}^{\prime} u_{j}^{\prime}}\right)=D_{i j}+\varphi_{i j}+G_{i j}-\varepsilon_{i j} $

(4)

式中：等号左侧第1项为应力时间项, 第2项为应力对流项, 右侧D_ij为湍流扩散项, φ_ij为压力应变项, G_ij为应力产生项, ε_ij为耗散项。

经过模化后的雷诺应力模型中的雷诺应力为二阶张量，需要求解6个雷诺应力分量方程和1个湍流耗散输运方程。RSM模型在求解雷诺应力时需要建立更多的方程，导致其计算效率低于2方程模型，但由于RSM模型直接求解雷诺应力张量，可以更加准确的反应雷诺应力的各向异性，一般来说在复杂湍流应用中RSM模型的适用性和精确性优于两方程模型^[15]。

由于在RPV整个流域中，平均雷诺数达到了10⁴量级，所以本文的多湍流模型方案选取SKE、RKE以及RSM这3个高雷诺数模型进行研究。上述3种模型有相同的壁面处理方式，可使用同样类型与密度的网格对多湍流模型混合计算精度与效率进行对比分析。MRANS计算时需要在区域边界面上对求解动量方程与湍流方程的重要参数数据进行传递，表 2显示了分段界面上各模型计算所需的参数。

图 9、10为应用不同的湍流模型的计算结果与实验值^[11]的对比，分别对应着DC不同高度(upper和lower)测点的混合标量值。

如图 9、10所示，在下降段上部2方程模型与RSM模型的计算结果差异很小且均与实验值符合较好。在下降段下部，2方程模型预测的标量场空间分布的坡度比实验值更加陡峭，而RSM明显与实验值的分布更加接近。可以看出，2方程模型明显的低估了下降段湍流交混的能力，这可能是将各向同性的2方程湍流模型应用在具有高的各向异性的结构上所导致的。同时，这种交混能力的低估在很多的有关FLOMIX的结果中被发现^[10]。

根据下降段的计算结果对比，本文基于SKE，RKE，RSM 3个模型进行互相匹配，设计多湍流模型耦合计算方案。因为多湍流模型计算方案的意义在于解决模拟大尺度部件时计算时间与计算精度难以平衡的问题，本文将遵循保证计算精度同时减小计算所需的资源的原则进行MRANS方案的设计。图 11为本文设计的MRANS计算方案，其中垂线连接的模型属于单一模型计算方案，斜线连接属于MRANS计算方案。通过各计算方案对整个流域的求解计算，可以从计算精度与效率两个方面综合评价MRANS方案与传统方案的区别。从表 2中可以看出，7方程模型包含有2方程模型所需要的全部参数，因此，在将7方程模型应用在2方程模型的上游时，无需补充参数。

3.2 MRANS方案验证

根据图 11的MRANS计算方案，分别对5种SKE, RKE, RSM的单一/混合模型进行计算。作为下腔室的入口，下降段出口处的流场及标量场结果的准确性会对下腔室流场的计算造成极大的影响。图 12分别展示了MRANS方案和单一湍流模型方案在下降段下部测点预测的速度和标量。图中黑色实线为实验值，虚线为实验值±30%范围。图 12(a)中速度为z轴方向速度。如图 12所示，由于MRANS模型在下降段使用RSM模型进行计算，因此在下降段下部测点的速度分布与RSM模型计算结果高度一致。SKE模型和RKE模型的计算结果也表现出一致性。在速度分布的计算中，RSM模型和MRANS模型在大部分区域与2方程模型精度相差不大，但在峰/谷值处与实验值更一致。

考虑到下降段标量场的影响范围，图 12(b)将标量未扩散的区域的数据予以剔除。相较于7方程模型，单一的两方程模型低估了低浓度区域的混合标量浓度，并且高估了高浓度区域的混和标量浓度。这与图 10的计算结果相符合。RSM的混合标量预测结果与实验值符合最好，而MRANS模型也仅在极少的位置与RSM模型有差异。无论是速度场还是标量场，MRANS模型都展现出了比SKE或RKE更优秀的精度。

图 13、14展示了不同RANS方案的RPV堆芯入口第16号和第25号通道的标量浓度的计算结果。在靠近下降段的第16号通道，单一的SKE和RKE模型预测的混合标量浓度上升非常迅速，在10 s内达到了充分混合的平稳阶段(在实验中上升过程约为15 s)，且过高的预测了平稳阶段的混合标量浓度。相较于2方程模型，RSM模型在上升阶段和平稳阶段的预测结果均与实验符合较好，两种MRANS方案中RSM-SKE模型在16通道均与RSM模型符合较好，且能较准确的描述标量随时间的变化过程。单一的SKE和RKE模型计算结果相类似，在上升阶段和平稳阶段预测结果与实验值的平均偏差分别为32%和17%，MRANS模型为11%和9%，RSM模型为5%和7%。而在相对更靠近中心的25号通道，单一的两方程模型预测的上升过程与实验符合较好，但仍低估了平稳阶段的标量浓度。所有的模型在上升阶段的预测与实验的平均偏差都在10%左右。在平稳阶段，单一的SKE和RKE模型的平均偏差为19%，MRANS模型为11%，RSM模型为5%。两种MRANS方案在上升阶段和平稳阶段的前半段预测符合较好，但在平稳阶段的后半段的混合标量浓度预测较RSM仍有较大偏差。之所以产生这种情况是因为，16通道为边缘通道，靠近下降段，流体从下降段进入下腔室后迅速进入通道，因此MRANS方案对该区域的准确性的提升较大。而25通道为中心通道，流体进入下腔室后需要流经分配孔板等几何结构后才能进入通道，因此MRANS方案对该区域的提升较小。

通过上述分析可以看出，基于RSM模型的MRANS方案可以得到与单一RSM模型符合较好的计算结果。由于不是全域使用更耗时的RSM模型，采用MRANS方案减少了大量的计算时间。表 3为本文所设计的MRANS方案与单一的RSM模型在Inter(R) Xeon(R) CPU E3-1280 v5，3.7 GHz，64 G的单机8核处理器条件下的计算效率。如表 3所示，MRANS方案相较于单一的RSM模型更加高效，在本文的网格规模下可以缩短约30%的CPU时间，且随着网格数量增加或LP网格数量在整体RPV模型的比重增加，可以增大这种MRANS方案的优势。

4 结论

1) 本文将棒束CFD分析中提出的分段方案运用到压力容器CFD分析中，设计了基于区域重叠方法的双向数据传递分段计算方案。

2) 分析了不同的湍流模型在各区域的计算结果，为分区计算选择合适的湍流模型。

3) 通过运用MRANS方案和所设计的双向数据传递分段计算方案构建了SKE-RSM与RKE-RSM两种MRANS方案，相比于传统的单一RSM湍流模型计算方案，本文提出的MRANS方案具有更高的效率。研究成果对提高RPV CFD仿真效率提供了一种新思路和途径。