自然对流是在热能转换与传递过程中一种常见的对流换热方式，由流体中的浮力引起，没有受迫速度^[1]。浮力是与密度成比例的物体力作用在有密度梯度的流体上产生的。流体的密度与温度有关，通常会随着温度的升高而降低，当重力作用在有温度梯度的流体上时，就会发生自然对流换热^[2-3]。自然对流独特的物理起因和换热特性，加上其流动过程中不需要外力做功，相比于强迫对流更具有经济性，使自然对流换热长期以来一直是传热学研究的重点、热点之一。对流换热通常需要解决的问题是确定对流换热系数，在自然对流中Nu的关系式可以直接计算对流换热系数，故而本文的核心点之一就是得到准确Nu的关系式。

在自然对流换热的研究中，广延空间内绕水平圆管的换热研究是传热学经典问题之一，其在换热器、核反应堆、太阳能收集器、冷凝器等工业设备里广泛存在。在传热学相关的教科书和期刊文献中，研究内容主要集中在单圆管和规则排列有间距多圆管的自然对流换热^[1-13]。其中，Morgan^[4]、Churchill和Chu^[5]在单圆管方面做了系统性的调查研究，针对不同的Ra、Pr范围给出对应Nu的关系式，在这之后的工业应用中普遍证实了这两篇文献中关系式的准确性和精确性。在工业实际应用之中圆管一般都是纵向、水平或者倾斜多个圆管较规则排列的，在这种结构中，相邻圆管之间的距离对流动和换热就会产生一定的影响，其换热特性并不是多个单圆管的简单叠加，或者线性关系的某种函数^[6-7]。基于这样的背景，竖直排列多圆管的绕流换热研究成为自然对流领域的另一个研究热点。为了揭示其流动换热特性，大量的研究者在过去的几十年里对其进行了详尽的研究^[6-12]。这些实验^[6-11]和数值模拟^[12]的适用条件，包括工质、阵列中的圆管数量N、圆管之间的中心距S及Ra或格拉晓夫数Gr的变化范围等，都已被简单地总结在表 1中(其中D为圆管直径，ΔT_L/ΔT_U为温差比)。综上所述，单圆管自然对流换热的研究可以为深入地研究多圆管自然对流换热机理奠定基础。对于竖直排列的多圆管，沿着重力方向最下面的圆管(即上游圆管)换热特性与单圆管非常接近，而上面圆管(即下游圆管)的换热能力则要受到S、表面温度T的影响。上述研究将多圆管阵列中下游圆管Nu的这种特别变化归因于预热影响和速度影响。

上述研究中，缺乏圆管之间间隙为零这样一种极限情况，即无间隙紧密接触。现在的换热器、核反应堆堆芯、高射炮炮管、枪的弹管普遍存在多圆管紧密接触这种极限情况，如图 1所示。紧密接触多圆管的换热特性目前还缺少全面、详细的研究，没有一个适用范围较强的经验公式。基于圆管自然对流换热的研究现状，并考虑其对于工业和军事的重要性，有必要对紧密接触品字形三圆管的自然对流换热进行全面的研究，以便为相关工业和军事应用提供精确的经验关系式，为其设计和优化提供参考。该研究的难点在于广延空间所导致的计算区域非常大，所需网格数量过多给计算带来麻烦。为解决这个难点，在之前的研究中，本课题组采用Bejan等^[13]提出的新计算模型，研究紧密接触水平和竖直排列双圆管的自然对流换热情况^[14-15]。研究结果表明，此模型可以大幅减小计算所需区域的大小，节约计算成本，解决长期困扰自然对流领域研究的大空间、高计算量的问题。在该研究的基础上及受到紧密接触三圆管枪(图 1)的启发下，本文研究水平放置紧密接触品字形三圆管(圆心成正三角形)在空气中的自然对流换热情况，研究范围为10≤Ra≤10⁶。

1 数理模型及计算方法

广延空间中物体绕流自然对流的研究一般是把所研究的物体放在一个足够大的人为设定的有限空间里，空间的边界对于该结构的流动传热需不施加任何的影响。因此，考虑把紧密接触品字形三圆管放在一个足够大的空间里，所有圆管都假设是定温的，其表面温度为一定值(T_w)，浸没在静止的空气中，周围空气的温度也设为一定值(T_∞)，并且T_w>T_∞。由于T_w>T_∞，圆管附近的空气与圆管表面存在一定的温度梯度，圆管表面附近的空气会被圆管表面加热，由于气体的热膨胀，在圆管附近被加热的空气密度将会降低，导致热空气与远处的冷空气之间存在一个密度梯度，进而在空气中产生浮力，在浮力的推动下，空气围绕圆管流动发生对流换热。一般情况下，为了使人为所设定的有限空间边界不影响圆管的流动及热传递，计算区域往往是非常大的，相应的网格量就会很大。虽然一般圆管绕流问题都会假设流动Z方向上没有分量，将本来的三维问题简化为二维问题，但是计算量一般还是无法承受的。巨大的计算区域不仅会使计算时间成本巨大，而且也会限制研究所覆盖的Ra范围，使获得一个能在较大Ra范围内有意义的传热经验关系式变得非常困难。本文计算模型如图 2所示，采用的方法是Bejan等^[13]提出, 假设流动是对称的，且是二维稳态的，重力方向为Y轴的负方向。另外，流体的热物性参数都用膜温T_m(=(T_w+T_∞)/2)的值，且除Y方向动量方程中的密度项外，其余变量都假设与温度无关^[14-15]。Y方向动量方程中的密度项假设仅随温度的变化而变化，它们之间的关系可由Boussinesq近似假设表示：

$ {\rho _\infty }-\rho \approx \rho \beta \left( {T-{T_\infty }} \right). $

(1)

在密度变化只由温度变化引起情况下其定义是

$ \beta =-\frac{1}{\rho }{\left( {\frac{{\partial \rho }}{{\partial T}}} \right)_p}. $

(2)

式中：ρ是工质的密度，ρ_∞是无穷远处的密度；β是体积热膨胀系数。显然，这个假设仅当区域中的温差很小时最为有效，小温差可以显著降低结果由热物性和β随温度变化而引起的误差。因此，为确保结果的准确性，本文中所有算例的温差都保持在1 K以内。由于流动的温差极小，β可以按照理想气体方程，将其设定为β=1/T_m。

假设工质是不可压缩的，并且忽略能量方程中的黏性耗散和辐射。控制方程组的数学表达式^{[1-3, 13-15]}如下所示，

连续性方程：

$ \frac{{\partial u}}{{\partial x}} + \frac{{\partial v}}{{\partial y}} = 0. $

(3)

X轴的动量方程：

$ u\frac{{\partial u}}{{\partial x}} + v\frac{{\partial u}}{{\partial y}} =-\frac{1}{\rho }\frac{{\partial p}}{{\partial x}} + \vartheta {\nabla ^2}u, $

(4)

Y轴的动量方程：

$ u\frac{{\partial u}}{{\partial x}} + v\frac{{\partial v}}{{\partial y}} =-\frac{1}{\rho }\frac{{\partial p}}{{\partial y}} + \vartheta {\nabla ^2}v + g\beta \left( {T-{T_\infty }} \right). $

(5)

能量方程：

$ u\frac{{\partial T}}{{\partial x}} + v\frac{{\partial T}}{{\partial y}} = \alpha {\nabla ^2}T. $

(6)

上述的控制方程中：u、v分别是x，y方向的速度；p是相对于无穷远处标准大气压p_∞的表压力；${\nabla ^2} = {\partial ^2}/\partial {x^2} + {\partial ^2}/\partial {y^2}$；g、$\vartheta $、α是重力加速度、工质的运动黏度和热扩散系数。

由于流动是关于垂直中心线对称的，所以本文的计算区域只需选取一半(X>0)作为计算区域。其边界条件的具体设置及无量纲参数请参见文献[13-15]。

上述控制方程在ANSYS Fluent (version 15.0.0)里通过有限体积法被离散化。压力-速度耦合格式采用压力耦合方程组的半隐式方法SIMPLE联立求解，动量和能量方程中的对流项采用二阶迎风格式进行离散化。为了获得更加准确的结果，连续性、动量和能量方程的收敛标准都设定为10^-6。然而，自然对流的传热计算往往是极难收敛的，想完全达到10^-6的收敛标准其计算时间相当漫长。因此，本文将圆管表面的平均热流密度以及阻力系数作为收敛的判定标准进行输出监视，当这两个参数的值保留4位有效小数位数且不随迭代步数变化时，就认为区域内的流动及换热已经达到稳定。

在自然对流的数值研究中，网格是影响计算结果的一个重要因素，对其合理的划分可以得到准确的结果，减小网格数量，提高计算效率。图 3是品字形三圆管计算的网格示意图，将其划分为两个区域，在圆管附近的区域I内，网格采用加密三角形非结构网格，有利于精确捕获边界层内的速度和温度的变化。而在边界层以外区域II则采用较稀疏的四边形网格，这样在不影响准确性的前提下，可以减少网格数量，进而提高计算效率。需要说明的是：本文三圆管附近两个网格之间的距离δ与圆管的直径D的比值，即为δ/D≤0.002 5, 这样就能满足网格质量而且保证计算结果的准确性^[13-15]。其中，A、B、C、D、E、F、O等点是为了后文分析方便表述而添加的。

2 计算结果及分析 2.1 模型方法验证

综上所述，为了获取品字形三圆管自然对流换热的准确结果，需要对数值模拟的模型和方法进行可行性与可靠性的验证。本文采用单圆管计算的结果与已有的实验或数值模拟的结果进行对比验证, 图 4给出单圆管Nu与Morgan^[4]、Churchill和Chu^[5]的公式对比图。研究结果符合很好，误差在6%以内，其详细的验证见本课题组已经发表的文献[14-15]。以上充分验证了模型方法、区域大小、网格精度的精确性。

2.2 温度场和流线场

在分析流动和传热问题的时候，借助于温度场和流线场可以很简单明了地表述研究问题的关键所在。图 5是三圆管附近的温度场(左侧)和流线场(右侧)随Ra变化的图示，X、Y轴的数值是无量纲处理后的，图注中0~1是无量纲温度。总体而言，随着Ra的增大，换热方式由热传导逐渐向对流换热转变，对流流动逐渐增强，温度和速度边界层随之变薄。圆管下面的旋涡(A点附近)是由圆管表面上的滞止点和圆管之间的接触点引起的：在浮力的驱动下，流体沿着中心线向上流动，在前驻点处速度降为零，并伴随着压力升高。从这点开始，流体沿着圆管表面分流，其中一边的流体在流到接触点时改变方向，向上游流去，然而来流的作用使流体又不能向上游回流，在这种情况下，旋涡在接触点下方形成。而沿着圆管表面向另一方向流动的流体则继续保持附着在圆管的表面上。从前驻点开始，压力沿表面逐渐降低，流体在顺压力梯度的作用下逐渐加速。但是在D点附近，流动受到后方圆管的高温预热及阻挡影响，流动产生分离：一部分开始回流，导致此处产生回流涡；一部分沿着后方圆管表面继续流动，压力最终会降至最小值，在圆管的后部会产生逆压力梯度。此时流体开始减速，随着流体的减速，表面上的速度梯度最终会在圆管表面后部的某一位置变为零，这个点称为分离点。在该点，圆管表面的流体动量不足以克服逆压力梯度，边界层发生分离，并在下游形成尾流。在Ra=10⁵和Ra=10⁶时，在其后方产生分离涡，这与单圆管相同，但是在Ra较小的情形时，分离涡产生得不明显，这反映出流动的增强对分离点的影响，在相同Ra的情况下堆积三圆管的流动强度小于单圆管的流动强度。由于在三圆管的空隙交汇处，与外界通道并没有对流换热，该区域始终维持在高温情形下，无需展示其温度场和流线场。

2.3 局部努塞尔数

局部努塞尔数Nu_θ的分布反应了圆管表面上每一点处温度的变化，是衡量圆管表面流动和传热特性的一个重要参数。其数值由于圆管表面的法向局部温度梯度决定，并且沿圆管表面连续变化。而圆管表面的温度梯度则主要受以下两个因素的影响：气体与圆管表面之间的温差及圆管表面的速度梯度。由于假设流动和温度分布是关于中心线对称的，因此图 6仅展示右边圆管在10≤Ra≤10⁶范围内Nu_θ的分布。需要说明的是，这里的度数是计算模型的弧长所对应的度数，在图 3三圆管的网格示意图中，前半圆管的弧长对应的度数为300°，即$\widehat {ABCD}$；后圆管圆弧对应的度数为150°，即$\widehat {DEF}$, 故而弧长总共对应的度数为450°。总的来说，Nu_θ均随Ra的增加而增加，这是由于作用在流体上的浮力随Ra的增加而增大，进而引起流动速度的增强，使得热边界层变得更薄，边界层热传导的热阻降低。相反，在较小Ra的情况下，流动较弱，速度梯度较小，传热的主要模式是热传导，流动对传热速率的贡献很小，因此其Nu_θ表面的变化很小，如图 6所示。在0~300°(即$\widehat {ABCD}$)其Nu_θ的值自0开始逐渐增加，前期在A点处正是堆积三管的滞止区，由于回流涡的存在，导致此处换热的能力很弱，所以Nu_θ的数值接近于0。在140°左右的时候出现Nu_θ的最大值，然后随着流动而减小，其Nu_θ的最大值与Ra呈正相关。在300°附近(即D点)，由于后圆管的预热与阻挡影响，形成回流涡，此处的换热能力趋近于0。在300°~450°(即$\widehat {DEF}$)，Nu_θ先是增加到一个较大的值，随后减小，但并没有趋近于0，由于其自身结构特性，在分离点以后受到羽流的影响较单圆管的要小，故而在分离点的后方形成分离涡，其换热能力变弱。

2.4 平均努塞尔数

平均努塞尔数Nu是衡量一个几何结构整体平均传热速率的重要参数，在实际工程的计算中通常是不可或缺的，Nu可以直接反映出圆管的整体换热能力强弱。图 7显示10≤Ra≤10⁶范围内单圆管和品字形三圆管Nu的变化趋势。可以看出，由于Ra的增加使对流增强，Ra与Nu呈现正相关。为了能够更准确地预测10≤Ra≤10⁶范围内水平紧密接触品字形三圆管结构的整体传热速率，参考文献[1-13]可发现，对于圆管Nu的拟合关联式普遍采用下面的形式：

$ \overline {Nu} = \frac{{hD}}{k} = CR{a^n}. $

(7)

C、n是不同的系数，针对于本文的模拟结果，拟合得到Nu的关联式如下

$ \overline {Nu} = 0.372\;5R{a^{0.223\;8}}. $

(8)

此拟合的曲线被绘制在图 7中。公式的可决系数R²=0.999 3，表明拟合公式与数值模拟结果之间的拟合程度很好，因此可以确定所获取的经验公式可以为水平紧密接触品字形三圆管结构的自然对流换热提供工程计算所需要的精确预测，同时，其公式为无量纲的变化因素拟合得到的，普适性较强。

此外，在图 7中，还对比了水平紧密接触品字形三圆管与单圆管的经验公式。其中，单圆管选取具有代表性的Morgan^[4]、Churchill和Chu^[5]的经验公式也是本文的验证公式。从图中可以看出，在相同Ra时，三圆管的Nu要比单圆管小得多。这是因为三圆管在紧密接触情况下，后方圆管的预热影响到前方圆管的流动换热，同时，前方圆管也会预热后方的圆管，这就导致三圆管边界层彼此相互影响，换热能力削弱，所以其Nu要比单圆管的小。而且，这种差异随Ra的增大而显著，偏差整体在55%~70%之间，在Ra=10⁶时差异约为70%。由此可见在二者之间存在这么大差异的情况下，如果仅仅用传统单圆管的经验公式来估算紧密接触品字形三圆管的传热速率，显然是不能满足工业计算要求的，这也体现出本文研究的重要性和必要性所在。

3 结论

本文利用CFD(计算流体力学)数值模拟的方式，在Pr=0.707 0，10≤Ra≤10⁶的范围内，得到紧密接触品字形三圆管Nu的关联式，进而可以计算该结构的对流换热系数。结果表明：

1) 随着Ra的增大，流动增强，对流换热能力增加，品字形三圆管的A、D点附近接触处的回流区也逐渐减小。在Ra=10⁵和Ra=10⁶的时候，在其后方产生分离涡，与单圆管相同，但是在Ra小的时候，产生的分离涡并不显著，这反映出流动的能力对分离点的影响。

2) Nu_θ的峰值出现在140°附近，且数值大小随着Ra的增加而增大。由于回流涡的影响，A、D处换热能力趋近于0，分离点后分离涡处，由于受到羽流的影响换热能力也很弱。

3) 由拟合的关系式可知Nu与Ra呈现正相关性，因品字形紧密接触三圆管的换热受到相邻圆管的影响，Nu要小于单圆管的Nu，这种差异随Ra的增大而显著，数值偏差整体在55%~70%之间，在大Ra的时候差异约为70%。