1. 引言
基于扰流的各类对流强化传热结构广泛存在于先进能源转换系统中，传热的强化通常伴随流动损失扩大的代价，因此，在结构设计中对流动损失进行精确预测，将有助于提升系统能效的优化效率。注意到，大多数强化传热结构在拓扑上具有特定的规则性和周期性，这为建立数学模型用于分析预测流动损失提供了便捷路径。本文以电子元器件冷却、航空发动机涡轮叶片冷却等紧凑式传热装置中常见的错排短柱阵列 [1] [2] [3] 为例，提出一种可用于定量分析和预测周期性规则拓扑流动损失的完备数学模型。
短柱通常用于受限空间内的强化传热，具有小于4.0的高径比(H/D) [4]，其流动传热特性不同于一般的列管式换热器 [5]。列管换热器中，由于端壁的影响可以忽略，流动近似二维。相比之下，短柱阵列中的流场包络于柱面和端壁，具有较强的三维性，使高径比的影响无法忽略。此外，柱间距(包括展向柱间距Sy/D和流向柱间距Sx/D)也对流动损失产生影响 [6]。事实上，影响规则短柱阵列流动和强化传热性能的结构因素有且仅有高径比和上述两个方向的柱间距，以往已就此展开过大量研究 [7] - [17]。然而，多数研究将侧重点放在结构的强化传热特性上，鲜有针对结构流动损失特性的分析。由对流强化传热的一般机理(掺混促进或边界层扰动)可知，传热强化通常带来流动损失(即流阻)的增加，对整体能效产生负面影响。因此对流阻进行定量的分析预测显得尤为重要。
一般将流阻以阻力系数的形式表达 [11] [15] [17] [18] [19] [20]，即：
(1)
其中ΔP是阵列进出口间的总压损失，ρ是流体密度，Umax是展向相邻柱体间最窄截面的通流速度，N是流向的柱排总数。需要指出的是，阻力系数的定义可能采用不同的形式，但定义式核心部分ΔP/ρU2是不变的，尽管有时用阵列进口速度(Uin)代替Umax，分母上的常数有时由2变为0.5 [21] [22] [23]。一些研究也用Darcy阻力系数表征流阻 [12] [13] [14] [16] [24] [25] [26]，即
(2)
其中，其中的Dh和L分别是阵列水力直径和流向长度。虽然阻力系数有不同的定义方式，但都需要对总压变化ΔP进行预测以表征流动损失。这就使经验方程的建立具有灵活性。表1给出了此前研究所得阻力系数经验方程。可以发现，多数方程都将结构参数的影响排除在外，仅考虑流动参数即Re数的影响 [22] [26] [27]，意味着此类方程仅适用于某一特定结构。Damerow [28] 和Jacob [29] 分别在经验方程中加入了流向间距Sx和展向间距Sy的影响，但其通用性依然受到限制，原因是高径比H/D和另一方向柱间距的影响仍被排除在外。主要原因是，此类经验方程的建立多基于有限的实验数据，数据体系的不完备降低了方程的普适性和准确性。
针对上述问题，本文通过数值模拟方法获得完备的流动阻力性能数据，基于此建立了新的经验方程。在此基础上，可对包括所有结构参数和流动参数在内的全部影响因素进行定量分析。整套方法可推广于任意具有周期性规则拓扑的流动损失分析及预测，具备较强的普适性。
Table 1. Empirical equations for friction factor of staggered pin fin array
表1. 错排短柱阵列阻力系数经验方程
2. 结果与讨论
2.1. 参数定义
如上所述，错排短柱阵列的全部结构参数有高径比(H/D)、流向柱间距(Sx/D)和展向柱间距(Sy/D)，图1给出了结构参数定义，其中D为短柱的直径。所建数学模型认为阵列流动损失取决于三种结构参数和雷诺数。雷诺数的确定以柱直径D作为特征尺度，以展向相邻柱间最窄截面的通流速度作为特征流速，则指数型经验方程具有如下基本形式：
(3)
这样，问题变为如何确定系数k和指数α、β、γ、δ的具体值。基于常见的应用场合和操作工况，分别为三种结构参数和雷诺数设定等间隔的5层值，如表2所示。在所考虑的雷诺数范围内，阵列中处于完全湍流状态。
Table 2. Definitions of structural parameters
表2. 结构参数定义
2.2. 数值方法设定
以往经验方程缺乏通用性的原因是：由于实验测量昂贵耗时，所得用于回归的性能数据不完备。针对这一问题，采用计算流体力学(CFD)方法代替实验，鉴于数值模拟的低成本和便捷性，可在短期内获得大量的性能数据，解决回归数据完备性的问题。由表2，可对共计125种结构进行分析，每种结构计算5种雷诺数，一共进行625次计算，可为模型建立提供充足的数据支撑。
周期规则拓扑为数值模拟的简化提供了便利，能够大幅减少计算负荷，提高数据收集的效率。对于短柱阵列，由于展向流场的重复性，可仅提取一个纵列，两侧分别赋予平移周期边界条件。进一步，由于沿高度方向流场的对称性，可提取端壁至中截面间的一半区域，中截面处赋予对称边界条件。这样，计算域的宽度恰为展向柱间距，高度恰为柱高的一半。尽管每排短柱引起的局部流动损失量级相当，为了取得统计意义，沿流向设置7排短柱。此外，在阵列进口前和出口后分别设置长度为10D的延伸段，以在进口处形成充分发展流动和防止出口回流。典型计算域和边界条件设定如图2所示。
进口给定速度，速度大小根据连续性法则由展向最窄截面的最大过流速度算得：
(4)
其中Umax由雷诺数计算获得：
(5)
出口给定静压，大小为1 atm，全部壁面给定恒温、速度无滑移边界。
Figure 2. Typical computational domain of staggered pin fin array (H/D = 2.5, Sx/D = 3.0, Sy/D = 3.5)
图2. 错排柱肋阵列典型计算域(H/D = 2.5, Sx/D = 3.0, Sy/D = 3.5)
数值模拟在基于有限体积法的ANSYS FLUENT 19.1软件中进行，通过求解雷诺平均的Naviers-Stokes方程，获得短柱阵列中的三维流场。以理想气体作为工作介质，假定粘度遵循Sutherland定律。利用ICEM CFD软件划分全域结构化六面体计算网格，实现计算域的空间离散。近壁面网格加密并保证Y+值不大于1。图3所示为部分结构的局部网格。
Figure 3. Examples of local computational grids
图3. 局部计算网格示例
为了在计算效率和经济性之间达到平衡，首先进行网格独立性检验，为典型阵列结构(H/D = 2.5, Sx/D = 3.0, Sy/D = 3.5)建立三套不同总数的网格系统。图4所示为计算所得总压损失随网格总数的变化，显然，在网格总数超过486288后，计算达到网格独立，因此以其对应尺度为全部125种结构划分网格。根据结构参数的不同，网格总数在344,008至665,600之间变化。
在进行大规模数值计算之前，首先需作算法验证，保证获得准确可靠的结果。参考Axtmann [11] 针对(H/D = 2.0, Sx/D = 2.5, Sy/D = 2.5)的阻力系数实测结果，对比不同的湍流模型，如图5所示。定性上，各模型均得出阻力系数随Re数减小并趋缓的趋势，定量上，standard k-ω模型取得了与实验数据最一致的结果，因此，后续的变结构多工况数值计算均采用standard k-ω模型。
Figure 4. Calculated pressure loss against mesh size (H/D = 2.5, Sx/D = 3.0, Sy/D = 3.5)
图4. 计算总压损失随网格总数变化(H/D = 2.5, Sx/D = 3.0, Sy/D = 3.5)
Figure 5. Variation of friction factor with Re number: Comparison between test results and turbulence models
图5. 阻力系数随Re数的变化：不同湍流模型与实测结果对比
2.3. 流动损失规律
图6所示为不同结构参数短柱阵列阻力系数随Re数的变化情况，从中可对高径比和柱间距的影响作初步定性分析。在考察单一结构参数的影响时，使其余两种结构参数分别固定于研究范围内的最小值(图6(a)、图6(c)、图6(e))和最大值(图6(b)、图6(d)、图6(f))，以求不失一般性。从图6(a)、图6(b)可以看出，阻力系数随高径比的减小而增大，原因显然是端壁影响增强加剧了三维流动掺混。高径比的影响在其本身较大时有所减弱，同样是因为流场有向二维衍进的趋势。这意味着，不能简单通过增加高径比来减小流动损失。
注意到，阻力系数大体随流向柱间距的扩大而增大，但并不总是单调的。由图6(c)可以看出，结构紧凑条件下相邻柱间沿程损失的影响盖过了柱本身造成的局部损失，因而出现阻力系数随Sx/D增大不减反增的现象。从图6(d)表明，结构稀疏条件下，最大阻力系数并未在Sx/D极值处获得，而是出现在Sx/D = 1.0和Sx/D = 5.0之间某处，表明沿程损失和局部损失之间存在相互平衡关系。对比图6(c)和图6(d)可知，受阵列部件间流动交互作用减弱的影响，Sx/D对流动损失的影响随阵列紧凑度的减小而减轻。
阻力系数随展向柱间距的扩大呈减小的趋势，原因是较大的Sy/D值意味着较小的阵列阻塞率。图6(e)显示，展向柱间距在阵列紧凑度较高时影响较大，原因同样是流动交互作用的加剧。值得注意的是，最大阻力系数在Sy/D = 2.5而非Sy/D = 1.5处获得，由此推断，展向紧凑度高到一定程度时，相邻柱间大尺度流动受到抑制，减小了动能耗散，从而减小流动损失。
总体来看，三种结构参数中，展向柱间距对流动损失的影响最明显，其次是高径比，流向柱间距的影响较小。尽管可通过扩大展向柱间距减小阵列流动损失，但需要同步考虑相邻柱间流动尺度抑制与否的影响。需要指出的是，图6只能为定性分析提供参考，为方便定量分析，须基于大量的流阻数据和对应的结构参数建立经验方程。
2.4. 数学模型建立
对短柱阵列流动损失的全部影响因素，包括结构参数H/D, Sx/D和Sy/D，以及流动参数雷诺数ReD，因此总压损失势必为上述因素的函数，即
(6)
若考虑单一因素对总压损失的影响，可以写出如下关系：
(7)
将方程组(7)可改写为
(8)
方程组(8)中各式可以合并为单一方程，即
(9)
方程(9)为多变量函数关系的通用形式，指明了基准函数值和任意目标函数值间的相互关系。一旦G函数的具体形式被确定，通过累乘即可由基准总压损失算得目标总压损失。
理论上，均匀阵列各排短柱造成的流动损失相等，阵列整体总压损失可用单排柱损失表征，即
(10)
将由数值模拟获得的总压损失数据拟合为指数形式，形成单一因素独立造成总压损失的影响，获得如下方程组
(11)
至此，每一影响因素基准值的选择将最终影响所形成经验方程的准确性和适用面。此例中，将各影响参数，取其考察范围内的中间值，即H/D = 2.5，Sx/D = 3.0，Sy/D = 3.5，ReD = 35,000。然后，在确定每一结构参数对应的指数时，其余两种结构参数分别取使阵列获得最高和最低紧凑度的水平。例如，确定高径比H/D对应指数α时，以Sx/D = 1.0和Sy/D = 1.5组合，形成最高的阵列紧凑度；以Sx/D = 5.0和Sy/D = 5.5组合，形成最低的阵列紧凑度。类似地，在确定雷诺数对应指数δ时，选取三种结构参数组合同样遵循使阵列获得最高(H/D = 0.5, Sx/D = 1.0, Sy/D = 1.5)和最低(H/D = 4.5, Sx/D = 5.0, Sy/D = 5.5)紧凑度的原则。如此设置可使经验方程的参数覆盖面拓宽。指数拟合结果为
(12)
式(12)中指数的绝对值可表征对应因素对流动损失的影响程度。能够看出，结构参数中，展向柱间距的影响程度最大，其次是高径比，而流向柱间距的影响最弱，这与图6的定性分析结论是一致的，也与Lawson等人 [9] 基于实测结果的分析一致。定量上，Sy/D的影响程度是Sx/D的12倍，是H/D的6倍。因此，调整展向柱间距是减小流动损失的最有效途径。
将式(12)代入方程(11)，并与方程(9)联立，可得：
(13)
方程(13)表示了任意参数组合和基准参数组合之间的一种恒等关系。基准参数选取中间值(H/D = 2.5, Sx/D = 3.0, Sy/D = 3.5, ReD = 35000)，可算得方程(13)右侧大小为1.2178 × 10−6，则排数为N的短柱阵列总压损失经验方程为
(14)
进一步，若以阻力系数表征，则表示为
(15)
值得一提的是，理论上，经验方程(14)和(15)的适用范围仍受所选参数范围的限制，准确性也受参数点密度的影响。但本文提出的是一种可行的方法，在解决具体工程问题时，可根据需要扩大参数范围，增加参数点密度，快速解决适用性和准确性的问题。本方法适用于绝大多数具有周期性规则拓扑强化结构的流动损失分析预测，本例中，影响因素包括三种结构参数和一种流动参数，但对于数学模型本身，参数的个数是灵活可变的。实际使用时，首先进行拓扑辨识，确定流动损失的影响因素，然后可根据大量的实测或数值模拟结果完成流动损失的定量分析预测。
3. 结论
本文提出了一种可用于定量分析预测周期性规则拓扑流动损失的数学模型，以典型错排短柱阵列为例，对模型的建立和应用进行了阐述，主要有如下结论：
1) 对于周期性规则拓扑，可充分利用结构特点进行数值简化，同时采用平移周期和对称面边界条件，快速获得大量性能参数的完备基础数据，为经验方程的建立和定量分析预测提供支撑；
2) 所建数学模型的基础数据可完全来自数值模拟，但模拟方法须经实测数据对比验证，以此确保数值计算结果的可靠性，进而可保证数学模型的准确性。正是由于数值模拟具有低成本的优势，所建数学模型具有基础数据更加完备的特点，可将全部影响因素纳入考虑范围，通用性远强于现存模型；
3) 通过所建立经验方程中指数的绝对值，可在定量层面上获得相应参数对流动损失的相对影响程度，为流动减阻设计和结构优化提供明确指导。
4) 所建数学模型的分析预测对象不限于流动阻力，而是适用于具有周期性统计意义多类性能指标，可将全部因素的影响用定量化的数学语言来表征，因而具有广泛的应用价值。
基金项目
江苏省高等学校自然科学研究面上项目(No. 19KJB470002)。
NOTES
*通讯作者。