, 卢文强 中国科学院大学工程科学学院, 北京 100049
2017年7月20日 收稿; 2017年10月20日 收修改稿
基金项目: 国家自然科学基金面上项目（51576189）资助
通信作者: 刘捷, E-mail:nauty@ucas.ac.cn
摘要: 采用Bejan提出的新计算模型，模拟计算水平堆积紧密接触品字形三圆管在空气中自然对流的换热过程。在普朗特数Pr=0.707 0，瑞利数10≤Ra≤106的范围内，得到该结构的稳态温度场和流线场、以及局部和平均努赛尔数随Ra变化的结果。在该结构中，两圆管之间相邻部分的空间会形成旋涡；局部努塞尔数Nuθ的最大值出现在140°附近，且随着Ra的增加而增大。同时，获得该结构的平均努塞尔数Nu随Ra变化的经验公式，进而可直接计算出自然对流换热系数，为工业应用提供参考。
关键词: 自然对流瑞利数努塞尔数数值模拟
Numerical simulation of laminar natural convection heat transfer from three horizontal attached cylinders
SHI Hongyan, LIU Jie
, LU Wenqiang College of Engineering Science, University of Chinese Academy of Sciences, Beijing 100049, China
Abstract: In this study, by applying the new computational model proposed by Bejan, laminar natural convection heat transfer from three horizontal attached cylinders has been studied numerically over the Rayleigh number range 10 ≤ Ra ≤ 106 and at the Prandtl number Pr=0.707 0. The comprehensive results consisting of streamlines, isothermal contours, and local and average Nusselt numbers are presented with different Rayleigh numbers. The obtained results show that a vortex is formed in the space between the adjacent two cylinders. The maximum local Nusselt number appears near θ=140°, and it increases with the Rayleigh number. Furthermore, this work accurately fits the equation dascribing the variation of average Nusselt number of the three cylinders with Rayleigh number. Hence the heat transfer coefficient of natural convection can be determined, and an accurate reference value is provided for industrial application.
Keywords: natural convectionRayleigh numberNusselt numbernumerical simulation
自然对流是在热能转换与传递过程中一种常见的对流换热方式，由流体中的浮力引起，没有受迫速度[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为圆管直径，ΔTL/ΔTU为温差比)。综上所述，单圆管自然对流换热的研究可以为深入地研究多圆管自然对流换热机理奠定基础。对于竖直排列的多圆管，沿着重力方向最下面的圆管(即上游圆管)换热特性与单圆管非常接近，而上面圆管(即下游圆管)的换热能力则要受到S、表面温度T的影响。上述研究将多圆管阵列中下游圆管Nu的这种特别变化归因于预热影响和速度影响。
Table 1
 表 1 竖直排列多圆管结构的研究现状Table 1 Research status of vertical arrangement multi-cylinders structure
| 表 1 竖直排列多圆管结构的研究现状Table 1 Research status of vertical arrangement multi-cylinders structure |
上述研究中，缺乏圆管之间间隙为零这样一种极限情况，即无间隙紧密接触。现在的换热器、核反应堆堆芯、高射炮炮管、枪的弹管普遍存在多圆管紧密接触这种极限情况，如图 1所示。紧密接触多圆管的换热特性目前还缺少全面、详细的研究，没有一个适用范围较强的经验公式。基于圆管自然对流换热的研究现状，并考虑其对于工业和军事的重要性，有必要对紧密接触品字形三圆管的自然对流换热进行全面的研究，以便为相关工业和军事应用提供精确的经验关系式，为其设计和优化提供参考。该研究的难点在于广延空间所导致的计算区域非常大，所需网格数量过多给计算带来麻烦。为解决这个难点，在之前的研究中，本课题组采用Bejan等[13]提出的新计算模型，研究紧密接触水平和竖直排列双圆管的自然对流换热情况[14-15]。研究结果表明，此模型可以大幅减小计算所需区域的大小，节约计算成本，解决长期困扰自然对流领域研究的大空间、高计算量的问题。在该研究的基础上及受到紧密接触三圆管枪(图 1)的启发下，本文研究水平放置紧密接触品字形三圆管(圆心成正三角形)在空气中的自然对流换热情况，研究范围为10≤Ra≤106。
Fig. 1
 | Download: JPG larger image |
| 图 1 紧密接触品字形三圆管 Fig. 1 The gun with three attached cylinders 图 1 紧密接触品字形三圆管 Fig. 1 The gun with three attached cylinders --> | |
1 数理模型及计算方法广延空间中物体绕流自然对流的研究一般是把所研究的物体放在一个足够大的人为设定的有限空间里，空间的边界对于该结构的流动传热需不施加任何的影响。因此，考虑把紧密接触品字形三圆管放在一个足够大的空间里，所有圆管都假设是定温的，其表面温度为一定值(Tw)，浸没在静止的空气中，周围空气的温度也设为一定值(T∞)，并且Tw>T∞。由于Tw>T∞，圆管附近的空气与圆管表面存在一定的温度梯度，圆管表面附近的空气会被圆管表面加热，由于气体的热膨胀，在圆管附近被加热的空气密度将会降低，导致热空气与远处的冷空气之间存在一个密度梯度，进而在空气中产生浮力，在浮力的推动下，空气围绕圆管流动发生对流换热。一般情况下，为了使人为所设定的有限空间边界不影响圆管的流动及热传递，计算区域往往是非常大的，相应的网格量就会很大。虽然一般圆管绕流问题都会假设流动Z方向上没有分量，将本来的三维问题简化为二维问题，但是计算量一般还是无法承受的。巨大的计算区域不仅会使计算时间成本巨大，而且也会限制研究所覆盖的Ra范围，使获得一个能在较大Ra范围内有意义的传热经验关系式变得非常困难。本文计算模型如图 2所示，采用的方法是Bejan等[13]提出, 假设流动是对称的，且是二维稳态的，重力方向为Y轴的负方向。另外，流体的热物性参数都用膜温Tm(=(Tw+T∞)/2)的值，且除Y方向动量方程中的密度项外，其余变量都假设与温度无关[14-15]。Y方向动量方程中的密度项假设仅随温度的变化而变化，它们之间的关系可由Boussinesq近似假设表示：
| ${\rho _\infty }-\rho \approx \rho \beta \left( {T-{T_\infty }} \right).$ | (1) |
 | Download: JPG larger image |
| 图 2 模拟模型示意图及边界条件 Fig. 2 Schematic representation of the flow and the computational domain 图 2 模拟模型示意图及边界条件 Fig. 2 Schematic representation of the flow and the computational domain --> | |
在密度变化只由温度变化引起情况下其定义是
| $\beta =-\frac{1}{\rho }{\left( {\frac{{\partial \rho }}{{\partial T}}} \right)_p}.$ | (2) |
假设工质是不可压缩的，并且忽略能量方程中的黏性耗散和辐射。控制方程组的数学表达式[1-3, 13-15]如下所示，
连续性方程：
| $\frac{{\partial u}}{{\partial x}} + \frac{{\partial v}}{{\partial y}} = 0.$ | (3) |
| $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) |
| $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) |
由于流动是关于垂直中心线对称的，所以本文的计算区域只需选取一半(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等点是为了后文分析方便表述而添加的。
Fig. 3
 | Download: JPG larger image |
| 图 3 计算区域的网格示意图 Fig. 3 Schematic representation of the computational grid 图 3 计算区域的网格示意图 Fig. 3 Schematic representation of the computational grid --> | |
2 计算结果及分析2.1 模型方法验证综上所述，为了获取品字形三圆管自然对流换热的准确结果，需要对数值模拟的模型和方法进行可行性与可靠性的验证。本文采用单圆管计算的结果与已有的实验或数值模拟的结果进行对比验证, 图 4给出单圆管Nu与Morgan[4]、Churchill和Chu[5]的公式对比图。研究结果符合很好，误差在6%以内，其详细的验证见本课题组已经发表的文献[14-15]。以上充分验证了模型方法、区域大小、网格精度的精确性。
Fig. 4
 | Download: JPG larger image |
| 图 4 平均努塞尔数结果对比 Fig. 4 Comparison of average Nusselt number values with the previous results 图 4 平均努塞尔数结果对比 Fig. 4 Comparison of average Nusselt number values with the previous results --> | |
2.2 温度场和流线场在分析流动和传热问题的时候，借助于温度场和流线场可以很简单明了地表述研究问题的关键所在。图 5是三圆管附近的温度场(左侧)和流线场(右侧)随Ra变化的图示，X、Y轴的数值是无量纲处理后的，图注中0~1是无量纲温度。总体而言，随着Ra的增大，换热方式由热传导逐渐向对流换热转变，对流流动逐渐增强，温度和速度边界层随之变薄。圆管下面的旋涡(A点附近)是由圆管表面上的滞止点和圆管之间的接触点引起的：在浮力的驱动下，流体沿着中心线向上流动，在前驻点处速度降为零，并伴随着压力升高。从这点开始，流体沿着圆管表面分流，其中一边的流体在流到接触点时改变方向，向上游流去，然而来流的作用使流体又不能向上游回流，在这种情况下，旋涡在接触点下方形成。而沿着圆管表面向另一方向流动的流体则继续保持附着在圆管的表面上。从前驻点开始，压力沿表面逐渐降低，流体在顺压力梯度的作用下逐渐加速。但是在D点附近，流动受到后方圆管的高温预热及阻挡影响，流动产生分离：一部分开始回流，导致此处产生回流涡；一部分沿着后方圆管表面继续流动，压力最终会降至最小值，在圆管的后部会产生逆压力梯度。此时流体开始减速，随着流体的减速，表面上的速度梯度最终会在圆管表面后部的某一位置变为零，这个点称为分离点。在该点，圆管表面的流体动量不足以克服逆压力梯度，边界层发生分离，并在下游形成尾流。在Ra=105和Ra=106时，在其后方产生分离涡，这与单圆管相同，但是在Ra较小的情形时，分离涡产生得不明显，这反映出流动的增强对分离点的影响，在相同Ra的情况下堆积三圆管的流动强度小于单圆管的流动强度。由于在三圆管的空隙交汇处，与外界通道并没有对流换热，该区域始终维持在高温情形下，无需展示其温度场和流线场。
Fig. 5
 | Download: JPG larger image |
| 图 5 在10≤Ra≤106圆管附近的温度场(左侧)和流线场(右侧) Fig. 5 Isotherm contours (left half) and streamline (right half) in the vicinity of the cylinders for 10≤Ra≤106 图 5 在10≤Ra≤106圆管附近的温度场(左侧)和流线场(右侧) Fig. 5 Isotherm contours (left half) and streamline (right half) in the vicinity of the cylinders for 10≤Ra≤106 --> | |
2.3 局部努塞尔数局部努塞尔数Nuθ的分布反应了圆管表面上每一点处温度的变化，是衡量圆管表面流动和传热特性的一个重要参数。其数值由于圆管表面的法向局部温度梯度决定，并且沿圆管表面连续变化。而圆管表面的温度梯度则主要受以下两个因素的影响：气体与圆管表面之间的温差及圆管表面的速度梯度。由于假设流动和温度分布是关于中心线对称的，因此图 6仅展示右边圆管在10≤Ra≤106范围内Nuθ的分布。需要说明的是，这里的度数是计算模型的弧长所对应的度数，在图 3三圆管的网格示意图中，前半圆管的弧长对应的度数为300°，即
Fig. 6
 | Download: JPG larger image |
| 图 6 在10≤Ra≤106时堆积三圆管表面Nuθ的分布 Fig. 6 Distribution of the local Nusselt number along the surface of the cylinders for 10≤Ra≤106 图 6 在10≤Ra≤106时堆积三圆管表面Nuθ的分布 Fig. 6 Distribution of the local Nusselt number along the surface of the cylinders for 10≤Ra≤106 --> | |
2.4 平均努塞尔数平均努塞尔数Nu是衡量一个几何结构整体平均传热速率的重要参数，在实际工程的计算中通常是不可或缺的，Nu可以直接反映出圆管的整体换热能力强弱。图 7显示10≤Ra≤106范围内单圆管和品字形三圆管Nu的变化趋势。可以看出，由于Ra的增加使对流增强，Ra与Nu呈现正相关。为了能够更准确地预测10≤Ra≤106范围内水平紧密接触品字形三圆管结构的整体传热速率，参考文献[1-13]可发现，对于圆管Nu的拟合关联式普遍采用下面的形式：
| $\overline {Nu} = \frac{{hD}}{k} = CR{a^n}.$ | (7) |
 | Download: JPG larger image |
| 图 7 模拟结果拟合的关系式和单圆管公式的比较 Fig. 7 The correlation for the present average Nusselt numbers and comparison of the obtained correlating equation with the previous correlations of single cylinder 图 7 模拟结果拟合的关系式和单圆管公式的比较 Fig. 7 The correlation for the present average Nusselt numbers and comparison of the obtained correlating equation with the previous correlations of single cylinder --> | |
C、n是不同的系数，针对于本文的模拟结果，拟合得到Nu的关联式如下
| $\overline {Nu} = 0.372\;5R{a^{0.223\;8}}.$ | (8) |
此外，在图 7中，还对比了水平紧密接触品字形三圆管与单圆管的经验公式。其中，单圆管选取具有代表性的Morgan[4]、Churchill和Chu[5]的经验公式也是本文的验证公式。从图中可以看出，在相同Ra时，三圆管的Nu要比单圆管小得多。这是因为三圆管在紧密接触情况下，后方圆管的预热影响到前方圆管的流动换热，同时，前方圆管也会预热后方的圆管，这就导致三圆管边界层彼此相互影响，换热能力削弱，所以其Nu要比单圆管的小。而且，这种差异随Ra的增大而显著，偏差整体在55%~70%之间，在Ra=106时差异约为70%。由此可见在二者之间存在这么大差异的情况下，如果仅仅用传统单圆管的经验公式来估算紧密接触品字形三圆管的传热速率，显然是不能满足工业计算要求的，这也体现出本文研究的重要性和必要性所在。
3 结论本文利用CFD(计算流体力学)数值模拟的方式，在Pr=0.707 0，10≤Ra≤106的范围内，得到紧密接触品字形三圆管Nu的关联式，进而可以计算该结构的对流换热系数。结果表明：
1) 随着Ra的增大，流动增强，对流换热能力增加，品字形三圆管的A、D点附近接触处的回流区也逐渐减小。在Ra=105和Ra=106的时候，在其后方产生分离涡，与单圆管相同，但是在Ra小的时候，产生的分离涡并不显著，这反映出流动的能力对分离点的影响。
2) Nuθ的峰值出现在140°附近，且数值大小随着Ra的增加而增大。由于回流涡的影响，A、D处换热能力趋近于0，分离点后分离涡处，由于受到羽流的影响换热能力也很弱。
3) 由拟合的关系式可知Nu与Ra呈现正相关性，因品字形紧密接触三圆管的换热受到相邻圆管的影响，Nu要小于单圆管的Nu，这种差异随Ra的增大而显著，数值偏差整体在55%~70%之间，在大Ra的时候差异约为70%。
参考文献
| [1] | Incropera F P, Dewitt D P, Bergman T L, et al. Fundamentals of heat and mass transfer[M]. 7th ed. Indiana: John Wiley & Sons, 2011: 595-630. |
| [2] | 杨世铭, 陶文铨. 传热学[M]. 4版. 北京: 高等教育出版社, 2006: 197-276. |
| [3] | Holman J P. Heat transfer[M]. 10th ed. New York: Mcgraw-Hill Education, 2010: 215-375. |
| [4] | Morgan V T. The overall convective heat transfer from smooth circular cylinders[J]. Advanced Heat Transfer, 1975, 11: 199-210. DOI:10.1016/S0065-2717(08)70075-3 |
| [5] | Churchill S W, Chu H S. Correlating equations for laminar and turbulent free convection from a horizontal cylinder[J]. International Journal of Heat and Mass Transfer, 1975, 18(9): 1049-1053. DOI:10.1016/0017-9310(75)90222-7 |
| [6] | Lieberman J, Gebhart B. Interactions in natural convection from an array of heated elements experimental[J]. International Journal of Heat and Mass Transfer, 1969, 12(11): 1385-1396. DOI:10.1016/0017-9310(69)90023-4 |
| [7] | Marsters G F. Arrays of heated horizontal cylinders in natural convection[J]. International Journal of Heat and Mass Transfer, 1972, 15(5): 921-933. DOI:10.1016/0017-9310(72)90231-1 |
| [8] | Sadeghipour M S, Asheghi M. Free convection heat transfer from arrays of vertically separated horizontal cylinders at low Rayleigh numbers[J]. International Journal of Heat and Mass Transfer, 1994, 37(1): 103-109. DOI:10.1016/0017-9310(94)90165-1 |
| [9] | Ashjaee M, Yousefi T. Experimental study of free convection heat transfer from horizontal isothermal cylinders arranged in vertical and inclined arrays[J]. Heat Transfer Engineering, 2007, 28(5): 460-471. DOI:10.1080/01457630601165822 |
| [10] | Grafsronningen S, Jensen A. Natural convection heat transfer from three vertically arranged horizontal cylinders with dissimilar separation distance at moderately high Rayleigh numbers[J]. International Journal of Heat and Mass Transfer, 2013, 57(2): 519-527. DOI:10.1016/j.ijheatmasstransfer.2012.10.076 |
| [11] | Kitamura K, Mitsuishi A, Suzuki T, et al. Fluid flow and heat transfer of natural convection induced around a vertical row of heated horizontal cylinders[J]. International Journal of Heat and Mass Transfer, 2016, 92: 414-429. DOI:10.1016/j.ijheatmasstransfer.2015.08.086 |
| [12] | Corcione M. Correlating equations for free convection heat transfer from horizontal isothermal cylinders set in a vertical array[J]. International Journal of Heat and Mass Transfer, 2005, 48(17): 3660-3673. DOI:10.1016/j.ijheatmasstransfer.2005.01.010 |
| [13] | Bejan A, Fowler A J, Stanescu G. The optimal spacing between horizontal cylinders in a fixed volume cooled by natural convection[J]. International Journal of Heat and Mass Transfer, 1995, 38(11): 2047-2055. DOI:10.1016/0017-9310(94)00312-J |
| [14] | Liu J, Liu H, Zhen Q, et al. Numerical investigation of the laminar natural convection heat transfer from two horizontally attached horizontal cylinders[J]. International Journal of Heat and Mass Transfer, 2017, 104: 517-532. DOI:10.1016/j.ijheatmasstransfer.2016.08.075 |
| [15] | Liu J, Liu H, Zhen Q, et al. Laminar natural convection heat transfer from a pair of attached horizontal cylinders set in a vertical array[J]. Applied Thermal Engineering, 2017, 115: 1004-1019. DOI:10.1016/j.applthermaleng.2017.01.029 |
