外延分区耦合瞬态计算流体力学计算技术

张历轩， 陈广亮， 田兆斐， 钱浩， 李瑞

(哈尔滨工程大学 核安全与仿真技术重点学科实验室，黑龙江 哈尔滨 150001)

堆芯的精细化瞬态计算流体力学(calculated fluid dynamics，CFD)是了解堆芯瞬态热工水力特性，发现堆芯薄弱环节，优化堆内构件设计的重要手段[1-2]。典型压水堆(pressurized water reactor，PWR)堆芯存在数万个狭窄流道，且堆芯内部包含千余层定位格架，每层格架含数百个刚凸、弹簧、搅混翼等具备弯折、扭曲特征的复杂结构[3]。针对压水堆中的5×5典型燃料组件，其需要0.9亿的网格量，400 G的内存资源。无法依靠常规的方法对其进行瞬态CFD分析[4]。

针对大规模堆芯稳态CFD仿真，Navarro等[5]在5×5棒束通道数值研究中，其将研究流域划分成7个子域，并使用压降作为参考验证分段耦合技术准确性，但未对子域的划分方法进行讨论。李小畅[6]针对轴向计算域的分段方式及子计算域出口边界条件类型对分段求解技术的数值模拟结果影响进行了讨论。陈广亮[4]针对5×5的棒束燃料组件通道，使用分段耦合方法完成了全高度的CFD稳态计算，并根据冷却剂通道不同区域的特征，采用不同湍流模型耦合来提高计算效率。然而，现阶段缺少分段耦合技术在瞬态CFD仿真中的应用研究。

本文基于稳态分段耦合技术，设计了瞬态分段耦合CFD计算方案，开发了瞬态调控程序，并通过外延分区优化了子域划分方案，使其适用于瞬态耦合，并使用阶跃型和周期性的入口条件对其进行验证，讨论了其适用性条件。

1 计算流体力学计算模型

1.1 物理模型及数值计算方法

选取压水堆燃料组件带搅混叶片、3跨高度的AFA-2G 1×2格架模型开展研究，格架模型如图1[4]所示。棒径为9.5 mm，栅距为12.6 mm，棒壁间距为2.55 mm,定位格架跨距520 mm。

图1 5×5 AFA2G 格架示意Fig.1 The structure for the 5×5 AFA 2G grid

对其划分子域的主要依据是不同堆内构件处的流场特性，光棒冷却剂通道(coolant subchannel, CS)是反应堆内冷却剂流动和换热的主要区域，流动阻力较小。定位格架(spacer grid, SG)是对燃料组件起到支撑和固定作用,其流动阻力较大。搅混翼(mixing vane, MV)则是坐落在部分定位格架上，它可以增加冷却剂的横向搅混[4]。

本文流动介质为水，入口温度为292.8 ℃，速度为3.76 m/s。设置周期性边界条件，壁面使用自定义函数加载热流密度，无滑移壁面，出口设置为自由流出。使用Fluent进行计算，压力-速度耦合选择SIMPLEC算法，湍流动能、耗散、能量项采用二阶迎风格式进行离散求解。能量方程收敛残差设置为1×10-6，其他项设置为1×10-4,采用标准壁面函数。

湍流模型选择雷诺应力模型(Reynolds stress model,RSM)。虽然RSM对计算资源需求较大[7]，但由于冷却剂在搅混翼区存在较强的横向搅混，具有较强的各项异性，采用RSM模型计算结果较为准确[8]。因此各子域间依次传递的信息不仅包括温度、压力、速度、湍动能、耗散率等，还应包括雷诺应力张量。RSM的控制方程推导如下所示：

湍流场中，各物理量均可以被分解成平均值与脉动值之和：

(1)

雷诺通过对纳维-斯托克斯方程(N-S)进行时间平均处理可得：

(2)

引入不封闭项-雷诺应力[9-10]。针对雷诺应力直接构造封闭方程，即为雷诺应力模型[11]。针对其所构造的输运方程为[12]:

(3)

针对右式各项分别建模，从而完成N-S方程的封闭。

将3跨高度1×2燃料组件从入口端至出口端，采用外延分区方案，按照光棒冷却剂通道、定位格架、搅混翼，依次划分为10个计算子域，具体子域划分方案见图2。针对3跨高度的1×2燃料组件和各个子域进行网格绘制。整体网格与各个子域网格均在光棒区使用六面体网格，在格架区使用四面体网格，设置3层边界层，y+值大约为100。

图2 外延分区子域划分方案Fig.2 The design of the simulation scheme using overlapping domain division

所有网格均使用压降进行网格无关性验证，如表1、2所示。对于整体计算网格，选择网格数量为175万的方案。ΔP为不同数量网格压降变化率，各子域也均采用压降进行网格无关性验证，以子域1为例，即冷却剂入口段，子域1选择网格数量为21万的方案。

表1 整体网格独立性验证Table 1 Overall grid independence verification

表2 子域1(CS)网格独立性验证

1.2 模型验证

为验证模型与计算方法的准确性，仿真结果与中国核动力设计研究院的实验数据[13-15]进行对比。选取距格架入口端50 mm和320 mm的下游处进行横流状态验证。在对堆芯进行CFD分析时，横流的状态[16]会影响冷却剂与壁面的换热以及流动阻力。常用的反映横流状态方法：1) 在关注截面上做一取样路径，对路径上检测点的横向速度进行速度分解；2) 将检测点的与取样路径垂直的横向速度与该截面的平均轴流速度做商；3) 以取样路径上各监测点的相对位置作为x坐标，以各点的商作为y坐标，绘制横向速度规律曲线。横向取样路径如图3所示。稳态工况下横流对比情况如图4所示，仿真结果与实验数据在趋势和数值上均比较接近，从而验证了模型的适用性。

图3 取样路径示意Fig.3 Sampling pathline for the analysis on lateral flow

图4 AFA2G格架轴向不同截面流速分布分析Fig.4 Analysis on the lateral flow status at different vertical cross-sections for AFA2G grid

1.3 瞬态分段耦合计算方案

稳态分段耦合方案是通过将研究区域划分为多个子域，依次传递边界条件从而实现全流域的CFD计算，其结果已得到验证[4-6]。

本文开发了调控程序，调用Fluent，使各子域在每个时间步上交换边界条件，顺次完成全流域的瞬态CFD计算。每个子域完成一个时间步长的计算后将出口参数传递给下一个子域，同时该子域继续进行瞬态计算，后续子域依次获取前一子域的边界条件，完成瞬态计算，直至整个流域瞬态仿真结束。

同时，该程序还具有保存各子域瞬态出口边界信息，调控特定子域参与计算和暂停等待的功能，既可以使全流域各个子域同步完成瞬态计算，也可以先完成前序部分子域的瞬态计算，保存各时间步的出口边界条件，以备后续子域调用，从而合理的规划CPU和内存的使用，实现多子域瞬态计算的分布并行。具体计算流程如图5所示。

图5 瞬态分段耦合计算方案Fig.5 Domain-divided solving technique in transient CFD

1.4 外延分区方案

在稳态分段耦合时增加了搅混翼格架下游计算域的长度保证了复杂流场计算结果的准确性[6]。然而将适用于稳态的子域划分方案应用至瞬态分段耦合中，其计算结果存在很大的误差。例如，在子域1的入口段引入流量阶跃后，监测子域2的出口截面平均流速，并将其与整体计算方案进行对比，如图6所示，其计算结果误差较大。

图6 原子域划分方案计算结果Fig.6 Results of the original strategy

因为相对于稳态计算，瞬态计算中搅混翼出口处不能保证为充分发展湍流，无法满足自由出口的应用条件[17-19]，即零扩散通量假设和质量平衡法则。因此在建模的过程中需要对其进行修正，将数值出口尽量设置到远离物理出口的位置，从而尽可能大的减小误差。

由此提出了外延分区的子域划分方案，重新进行了几何建模以及网格划分，将原子域的出口部分进行了延长，并重新绘制了网格。子域间数据传递时，将物理出口的边界条件传递给下一子域的进口。

尝试将原子域延长自身轴向长度的5%、10%、20%。对比其均方根误差和计算所需时间来选择子域延长方案。本文将整体计算方案仿真的数值作为真实值，并以此来计算外延分区耦合方案的误差率和均方根误差(root mean square error, RMSE)：

(4)

(5)

以格架搅混翼区为例，对比在周期性入口条件下，不同外延方案计算1 000个时间步出口截面轴流速度的RMSE和计算耗时。

经对比，延长子域轴向长度的20%可以降低与整体计算方案相较的RMSE，同时其计算耗时也不会大幅增加。因此本文的计算模型均采用外延子域自身20%轴向长度的计算方案。

表3 不同外延方案效率对比Table 3 Analysis of the performance of different strategies

2 瞬态计算适用性分析

2.1 流量阶跃响应特性分析

本文使用阶跃型和周期性入口边界条件对外延分区耦合瞬态CFD计算方案的有效性进行了验证。稳态时，将子域1的入口速度设置为3.76 m/s,所有子域都计算收敛后，将其入口的流量降低20%，调控所有子域进行瞬态计算。由图7可知，通过外延分区后，整体计算方案与分段耦合方案在瞬态仿真时，子域2出口的平均轴向流速的误差率在整个仿真时间内都很小，不到千分之一。相较于图6计算结果有了很大改善。

图7 外延分区耦合方案Fig.7 Results of the overlapping domain division

同样采用横向速度与轴向速度比值来分析横流特性，取样位置如图2所示。当引入流量下降的阶跃入口边界条件时，子域6(MV)出口截面的横向速度各时刻分布规律变化如图8所示，可以看到降低流量时，搅混翼出口横流速度与轴流速度的比值快速增大。2种方案均能反应这一特性，并且2种方案所反映的横流特性也基本一致，存在的误差从初始稳态至瞬态计算过程中未增加，从而验证了外延分区耦合瞬态CFD计算方案在阶跃型边界条件下仿真横流状态的准确性。

图8 不同时刻搅混翼出口横流特性分析Fig.8 Analysis on the lateral flow status at different times

2.2 横向边界适用性分析

在流域横向使用周期性条件对外延分区耦合方案进行测试，其初始稳态与流量阶跃案例一致。通过对子域1的入口流量引入周期为0.025 s的20%流量波动，在瞬态计算中，监测各个子域的出口截面平均流速，并与相同边界条件的整体计算方案进行对比，各横向截面对比结果如图9所示。各子域的划分方案如图5所示。

图9 不同子域出口轴向流速Fig.9 Axial velocity at the outlet of different subdomains

以整体计算方案参数为真实值进行对比，在子域1的出口，其平均流速的误差率仅有1×10-4，效果良好，如图9(a)所示。随着耦合子域数量的增加，后序子域的误差逐渐增大，最大误差率不超过7×10-3，最大RMSE不超过2×10-2。

图10比较了经过一个周期后，t=2×10-5s时刻外延分区耦合方案与整体计算方案的压降与温度在轴向上的分布，其温度分布在轴向上吻合良好，压降的数值和趋势也基本一致。

通过对比横向截面与轴向的参数，验证了外延分区耦合瞬态CFD计算技术在阶跃性和周期性边界条件下的适用性与准确性。在稳态分段耦合计算中，出口边界参数法向梯度值较小，以光棒区出口为例，其稳态计算中，如表4所示，出口截面参数的法向梯度相较于瞬态过程小很多。因此基于物理量的局部单相化假设条件[20]，以及出口零法向梯度的充分发展边界条件[21]，将出口设置为自由流出满足数值计算条件要求，计算结果准确。将出口的边界条件条件依次传递给下游子域，仍然能够保证计算结果的准确性。而瞬态变化过程中，热工参数波动非常剧烈的，出口边界参数法向梯度很大，因此各子域出口处不能保证为充分发展的湍流，无法满足自由流出条件要求，导致仿真结果失真。

图10 t=2×10-5 s时刻整体计算方案与外延分区耦合方案轴向压降与温度对比Fig.10 Comparison of axial parameters between overlapping domain division and overall calculation scheme at 2×10-5 seconds

表4 子域1出口处参数的法向梯度

而通过外延分区后，将数值出口设置在远离物理出口的地方，减少了不合理出口边界条件给计算引入的误差。

3 结论

1)本文研究了高适用且高效的瞬态CFD计算分析技术，提出的外延分区耦合瞬态CFD计算分析技术可以实现调控CFD程序实时切换边界条件，优化子域划分方案。

2)经对比分析验证该方案在瞬态计算中横截面速度、横流状态以及轴向温度、压降分布等均与整体计算方案吻合较好，具备工程适用性，并且降低了单次计算负荷，解决了堆芯大流域瞬态CFD计算对CPU和内存要求过高的问题。