迷宫型节流件压降的数值模拟研究

2020-06-21陆道纲秦亥琦钟达文

核科学与工程 2020年2期

宋怡，陆道纲，秦亥琦，钟达文，曹琼

(华北电力大学非能动核能安全技术重点实验室，北京 102206)

池式钠冷快堆主冷却系统的流道包括主循环流道和五个分支流道。在分支流道中，迷宫型节流件被安装在一些重要的部件前后来调节压降和控制流量，例如在主循环泵的旁通流量计中，钠从泵后侧引出的接管流出，依次流经迷宫型节流件、连接管道、流量计、出口管，最终流进主循环泵的泵箱。除此之外，泵支承冷却系统、堆容器冷却系统、堆内电离室冷却系统中也用到了迷宫型节流件。当冷却剂钠流过这些迷宫型节流件时，钠温均为360 ℃[1]。由于在设计和运行过程中反应堆的热工水力参数经常被调整，因此通过这些迷宫型节流件的流量和需要实现的压降也会有所调整。为了更加合理地挑选迷宫型节流件，需要较为准确地评估它的压降。

所有三维迷宫型节流件均由厚度相同、高度相同且高于管道半径的节流片组成，相邻节流片间距相同且呈180°交错排列。节流件放置在圆形管道中，示意图如图1所示。

图1 迷宫型节流件示意图

评估压降通常采用实验的方法，但是实验不仅费时费力，而且难以观察流动细节，近年来，数值模拟已经成为计算压降的常用方法。国内外许多学者用数值模拟方法研究了各种节流件的流动和阻力特性。Abdulrazaq A.使用realizablek-ε模型研究了两个相似的连续排列的斜边孔板的流动特性。发现收缩位置与第一个孔有关，与第二个孔无关，且随孔径比的变化而变化[2]。Ehsan S应用数值模拟方法，优化涡流管的参数和几何形状，获得了能够实现最高制冷效率的节流件锥长[3]。Guohui Gan采用CFD方法预测了方形薄板上单孔和多孔的压力损失特性，结论是多孔的压力损失系数比单孔高约0.5[4]。Ivan Okhotnikov对一种新型单级节流阀进行了CFD数值模拟，预测了由于稳定流动引起的扭矩和压降，分析了开口面积、流量系数和射流角度的图像，证明了该设计在高流量的流体动力系统中使用的有效性[5]。王海民通过数值模拟方法研究了由一系列直角转弯组成的曲折迷宫通道中的流速和压力分布，证明了由串联通道产生的多级压降可以在控制流速合理的同时防止空化现象的产生[6]。Han Xu模拟了带有环形塞的止回阀内的湍流流动，通过与实验结果对比，证明RNGk-ε模型足以预测该阀的水动力特性[7]。

尽管数值模拟很便捷，但它的准确性受很多因素影响，尤其是湍流模型的选择。目前最适合迷宫型节流件数值模拟的湍流模型尚不明确，数值模拟的结果也缺少实验验证。

另一方面，迷宫型节流件的压降是由多种因素决定的，包括相邻节流板的间距、节流板的高度、节流板的数目和流过节流件的流量。目前关于迷宫型节流件的压降和这些因素之间的关系尚不明确。多元回归分析方法是变量关系研究常用的方法，通过它拟合出显式公式不仅可以较为直观地表明变量之间的关系，而且可以通过已有数据来预测更多数据。Qi Sun使用计算流体力学方法，分析了颗粒雷诺数、流动角和阻力系数百分比对高温气冷堆中简化的石墨粉尘颗粒阻力系数的影响。并基于Rosendahl提出的框架采用多元回归方法拟合了更精确的阻力系数公式。他还应用最小二乘法求出相关参数，显著降低了预测阻力系数的误差[8]。Karim Nassar对海堤的水动力性能进行了实验，分析了各种因素对海堤水动力性能的影响，利用多元非线性回归方法开发了新的经验公式。不仅用于理解拟建海堤的水动力特性，并且预测了某些特性。基于实验结果，这些公式用软件SPSS进行了验证[9]。

基于上述认识，本文将研究计算迷宫型节流件压降的数值模拟方法，并且通过数值模拟结果与实验结果的对比，选择最可信的湍流模型。本文将使用经过验证的数值模拟方法，分析迷宫型节流件的压降随各种因素变化的趋势，最终通过多元回归方法拟合压降与这些因素之间关系的经验公式。

1 迷宫型节流件的数值模拟方法及其验证

1.1 迷宫型节流件的数值模拟方法

将通过实验验证的三维迷宫型节流件，放置在内径123 mm的圆形管道中。它由两块节流片组成，每块节流片厚度为7 mm，高度为78 mm，相邻节流片间距为16 mm，呈180°交错排列，如图2所示。为了使数值模拟的几何模型与实验的几何模型完全一致，我们将节流件的入口延长615 mm，出口延长1 230 mm。

迷宫型节流件的几何模型在Solid Works中创建，然后将其导入ICEM CFD生成非结构网格，并在近壁区生成良好的边界层，生成的网格将导入Fluent进行计算。Fluent的计算过程基于有限元方法，也就是将连续的流体域划分为离散单元，再对每个单元求解方程。如果单元数太少，难以模拟流动细节，计算结果将产生较大误差；如果单元数太多，计算机负担过重，会影响计算速度。因此，在选择湍流模型以前，我们分析了网格独立性，确定了能够保证计算精度的最小网格数。网格独立性分析(使用标准k-ε模型，雷诺数为852 744)的结果如图3所示，220万单元数的网格被选为最佳网格。

图2 数值模拟的几何模型

图3 网格独立性分析

由于目前最适合迷宫型节流件数值模拟的湍流模型尚不明确，因此我们运用了多个湍流模型来重复数值模拟过程，包括标准k-ε模型、标准k-ω模型、可实现的k-ε模型,RNGk-ε模型,SSTk-ω模型和Spalart-Allmaras模型。它们的控制方程可以从已经发表的文献中查阅[10-18]。将这些湍流模型计算得到的数值模拟结果与实验结果对比，得到最适合迷宫型节流件的湍流模型。

数值模拟中使用的流体是360 ℃的液态钠，与快堆主冷却系统支路中的钠温一致。液钠被认为是单相的、不可压缩的，它的物理性质不随温度变化。Fluent中的边界条件如表1所示。

表1 边界条件

1.2 验证

1.2.1 实验介绍

实验设备包括水箱、泵、阀门、加热器、流量计、差压计和节流件的实验段，如图4所示。其中水箱用来储存纯净水，并为实验环路提供水源。加热器用来维持环路中的水温。

图4 实验装置和环路示意图

节流件的实验段用不锈钢锻造，并通过法兰连接到实验环路。节流件的前后均布置有压力引管，差压计通过压力引管连接到实验环路。为了确保差压变送器所测的压降不受流量分布不均的影响，我们将取压点布置在相同高度，且采用环形四点取压方式。实验中，流体从节流件底部的入口流入，流经节流件后从顶部的出口流出，通过调节阀门的开度来改变流量。流量计和差压计自动收集体积流量和压降的数据，经信号转换模块转换后显示并储存在电脑端。

鉴于液态钠价格高昂及其危险性，我们使用水来代替钠进行实验。为了能够将实验结果应用于原始模型，实验需要满足相似理论，包括几何相似、运动相似和动力相似。

采用84 ℃的水进行实验，得到压降随体积流量变化的实验结果如图5所示。

图5 压降随体积流量变化的实验结果

1.2.2 数值模拟结果与实验结果的对比

我们使用不同湍流模型重复进行了迷宫型节流件的数值模拟过程，质量流量范围是11～31 kg/s，雷诺数范围是453 235～1 199 368。数值模拟结果与实验结果的对比如图6所示。在整个流量范围内，使用标准k-ε模型得到的数值模拟结果与实验结果的误差均在3%以内。因此，我们选择标准k-ε模型为计算迷宫型节流件压降最可信的湍流模型。

2 迷宫型节流件水力特性的分析

使用上述经过验证的数值模拟方法，我们分析了迷宫型节流件的水力特性，节流件由两到四块节流片组成，放置在内径为100 mm的管道中。

图6 数值模拟与实验得到的压降对比

2.1 空化现象及发生条件

在钠冷快堆中，一旦在任何部件内出现气泡，气泡就会混杂在液钠中，随液钠流向其他部件，这会恶化冷却剂的传热性能，对反应堆造成潜在的危害。当液钠的压力低于该温度下的饱和蒸气压时，就会出现空化现象并产生气泡。在节流件中，由于局部压力的快速下降，空化现象更容易发生。由于迷宫型节流件被安装在钠冷快堆的主冷却系统中，因此流过它的液钠的绝对压力应当高于360 ℃下的饱和蒸气压，防止空化现象的发生。

液态钠的饱和蒸气压与温度的关系可以通过多种方法获得。Golden和Tokar对这些方法进行了分析讨论[19]并获得了温度低于1 250 K时的最佳方程，该方程由Ditchburn和Gilmour给出[20]：

(1)

式中：Ps——饱和蒸气压,Pa;

T——温度,K。

当温度为360 ℃时，根据式(1)可以计算出饱和蒸气压为1.36×10-4Pa。

在迷宫型节流件中，流体每流过一个节流片，压力都会迅速下降然后缓慢回升。另外，除了最后一个节流片以外，每个节流片的后侧都有一个小漩涡，漩涡中心的压力低于周围区域的压力。因此压力最小的位置位于倒数第二块节流板的后侧。(见图7的A、B、C点)

图7 流线图与绝对压力云图

2.2 压降随各因素变化的趋势

迷宫型节流件的压降是由很多因素决定的，包括相邻节流板的间距(D)、节流板的高度(H)、节流板的数目(N)和流过节流件的流量(M)，如图8所示。本文分别分析压降随这些因素变化的趋势。

图8 影响迷宫型节流件压降的各因素

在迷宫型节流件中，H和D决定了流体的流通面积，因此我们首先分析压降与H、D之间的关系。随着H减小或D增大，流通面积都会增大，进而压降减小。当流体面积大到一定程度再继续增大时，压降几乎不变。另外，当流通面积太小时，空化现象更易发生。

图9展示了压降随D、H的变化趋势。当D小于引起空化的临界间距时，液钠的最小绝对压力会低于该温度下的饱和蒸气压，将导致空化并造成危险；当D小于25 mm，随着D增加，压降显著减小；当D大于25 mm时，随着D增加，压降几乎不变。因此，D的最佳范围是从引起空化的临界间距(Dmin)到25 mm。当H小于65 mm时，随着H增大，压降缓慢增加；当H大于65 mm时，随着H增大，压降显著增加。

图9 压降随D、H的变化趋势

图10展示了压降随N、M的变化趋势。由于迷宫型节流件的降压形式是多级降压，因此随着N增加，压降逐渐增大。随着M增加，单位时间内流过节流件的流体体积增大，因此压降增大。

图10 压降随N、M的变化趋势

2.3 迷宫型节流件压降的经验公式

由于在设计和运行过程中反应堆的热工水力参数经常被调整，因此通过这些迷宫型节流件的流量和需要实现的压降也会有所调整。我们定量分析了大量数值模拟结果，用多元非线性回归分析方法拟合了经验公式，拟合出的经验公式可以被用作评估迷宫型节流件压降的参考。经验公式的适用范围是管道内径100 mm，H为60～80 mm，D为Dmin—25 mm，N为2～4，M为4～12 kg/s。

由于四个自变量(N,M,H,D)在物理意义上相互独立，因此我们分别画出了表示压降P与这些自变量之间关系的散点图，如图11所示。

图11 压降与自变量关系的散点图

如图11所示，压降与N、M呈现一定的线性关系；此外，由于N和M维度的样本数量有限，限制了公式在这些维度上的拟合能力。因此，我们假定压降与N、M的关系是线性的。

如图11所示，压降与H和D的关系呈现非线性的高幂次关系；如图9所示，当其他自变量被固定，压降随H和D的变化趋势都大致呈现二次关系。因此，我们假定压降与H和D的关系是二次的。图12展示了不同M下压降随H和D变化的散点图，进一步说明了我们假设的二次关系是合理的。

图12 不同M下压降随H、D变化的散点图

根据上述分析，我们提出拟合函数的形式为：

P=a×N+b×M+c2×H2+c1

×H+d2×D2+d1×D+e

(2)

式中a,b,c2,c1,d2,d1,e均为待定系数。

对参数的求解属于回归分析中的参数估计问题，通常采用最小二乘估计方法，最小二乘法是一种通过最小化误差平方和寻找最佳函数匹配的优化算法，用公式化表示为：

(3)

其中ym是拟合曲线上的点，yi是观测即实验得到的真实值，对于一元线性回归y=a+bx，为使上式最小可分别对a、b求偏导，令偏导等于0可计算出a、b的值。可以通过数据矩阵、参数矩阵和观测值矩阵获得参数的显式解，计算过程如下：

将xi,yi(i=1,2,…,n)代入y=a+bx中：

(4)

令

本文中的多元非线性回归可以转换成线性回归方程求解。使用同样的方法可以求出系数矩阵：

(5)

我们基于计算机机器学习工具库sk-learn编程辅助计算，估计出各项系数的可靠经验值。我们使用了472个数值模拟得到的数据并随机分为训练集和测试集，其中测试集的数据占数据总数的20%。训练集的作用是计算拟合公式的系数，测试集的作用是对训练集拟合出的显式函数结果进行评估。通过计算，我们得到的经验公式为：

(6)

式中：P——迷宫型节流件的压降,kPa;

N——节流片的数目;

M——流过节流件的质量流量,kg/s;

H——节流片的高度,mm;

D——相邻节流片的间距,mm。

考虑到压降的物理意义(不会出现负值)，我们对公式(3)进行了规范，得到新的函数表达式为：

P′=Relu(P)

(7)

式中：P′——修正过的压降;

Relu——线性整流函数，它的作用是把负数都变为0而正数保持不变。

为了对拟合出的公式精度即回归分析的效果进行评估，我们选择了统计学中对非线性拟合常用的评估指标决定系数R2，其含义为由自变量变化引起的因变量变化占总体因变量变化的比例，该值越大证明函数拟合性能越优。其计算方法为：

平均观察值:

(8)

于是可以得到总平方和:

(9)

回归平方和和残差平方和分别为:

(10)

(11)

最终决定系数被定义为:

(12)

我们用测试集对训练集拟合得到的函数进行统计学评估，计算得到决定系数R2为0.888。在一般的统计学回归分析应用问题中，通常认为R2大于0.8即可认定函数拟合精度高，因此我们认为通过回归分析计算得到的经验函数合理有效[21]。同时，拟合出的经验公式通过了关于各自变量的显著性检验，说明了各自变量都对压降有着重要的影响。此外，该经验公式对数值模拟数据的均方根误差为32.259，比对其他函数形式模型(包括N、M、H、D的高次幂模型和H、D的低次幂模型)的均方根误差都要小，说明函数既没有欠拟合也没有过拟合。我们画出了测试集真实数据和公式输出数据的对比图像，如图13所示，可以较为直观地看出拟合到的函数是可信的。